pktools  2.6.7
Processing Kernel for geospatial data
pkcomposite.cc
1 /**********************************************************************
2 pkcomposite.cc: program to mosaic and composite geo-referenced images
3 Copyright (C) 2008-2014 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include <algorithm>
21 #include <vector>
22 #include <iostream>
23 #include <string>
24 #include "imageclasses/ImgReaderGdal.h"
25 #include "imageclasses/ImgWriterGdal.h"
26 #include "imageclasses/ImgReaderOgr.h"
27 #include "base/Vector2d.h"
28 #include "base/Optionpk.h"
29 #include "algorithms/StatFactory.h"
30 
31 /******************************************************************************/
138 namespace crule{
139  enum CRULE_TYPE {overwrite=0, maxndvi=1, maxband=2, minband=3, validband=4, mean=5, mode=6, median=7,sum=8,minallbands=9,maxallbands=10,stdev=11};
140 }
141 
142 using namespace std;
143 
144 int main(int argc, char *argv[])
145 {
146  Optionpk<string> input_opt("i", "input", "Input image file(s). If input contains multiple images, a multi-band output is created");
147  Optionpk<string> output_opt("o", "output", "Output image file");
148  Optionpk<int> band_opt("b", "band", "band index(es) to crop (leave empty if all bands must be retained)");
149  Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (empty: keep original resolution)");
150  Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (empty: keep original resolution)");
151  Optionpk<string> extent_opt("e", "extent", "get boundary from extent from polygons in vector file");
152  Optionpk<bool> cut_opt("cut", "crop_to_cutline", "Crop the extent of the target dataset to the extent of the cutline.",false);
153  Optionpk<string> eoption_opt("eo","eo", "special extent options controlling rasterization: ATTRIBUTE|CHUNKYSIZE|ALL_TOUCHED|BURN_VALUE_FROM|MERGE_ALG, e.g., -eo ATTRIBUTE=fieldname");
154  Optionpk<string> mask_opt("m", "mask", "Use the first band of the specified file as a validity mask (0 is nodata).");
155  Optionpk<float> msknodata_opt("msknodata", "msknodata", "Mask value not to consider for composite.", 0);
156  Optionpk<short> mskband_opt("mskband", "mskband", "Mask band to read (0 indexed)", 0);
157  Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box", 0.0);
158  Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box", 0.0);
159  Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box", 0.0);
160  Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box", 0.0);
161  Optionpk<string> crule_opt("cr", "crule", "Composite rule (overwrite, maxndvi, maxband, minband, mean, mode (only for byte images), median, sum, maxallbands, minallbands, stdev", "overwrite");
162  Optionpk<int> ruleBand_opt("cb", "cband", "band index used for the composite rule (e.g., for ndvi, use --cband=0 --cband=1 with 0 and 1 indices for red and nir band respectively", 0);
163  Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "invalid value(s) for input raster dataset");
164  Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Band(s) in input image to check if pixel is valid (used for srcnodata, min and max options)", 0);
165  Optionpk<double> minValue_opt("min", "min", "flag values smaller or equal to this value as invalid.");
166  Optionpk<double> maxValue_opt("max", "max", "flag values larger or equal to this value as invalid.");
167  Optionpk<double> dstnodata_opt("dstnodata", "dstnodata", "nodata value to put in output raster dataset if not valid or out of bounds.", 0);
168  Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation).", "near");
169  Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
170  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
171  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
172  Optionpk<string> projection_opt("a_srs", "a_srs", "Override the spatial reference for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid");
173  Optionpk<short> file_opt("file", "file", "write number of observations (1) or sequence nr of selected file (2) for each pixels as additional layer in composite", 0);
174  Optionpk<short> weight_opt("w", "weight", "Weights (type: short) for the composite, use one weight for each input file in same order as input files are provided). Use value 1 for equal weights.", 1);
175  Optionpk<short> class_opt("c", "class", "classes for multi-band output image: each band represents the number of observations for one specific class. Use value 0 for no multi-band output image.", 0);
176  Optionpk<string> colorTable_opt("ct", "ct", "color table file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
177  Optionpk<string> description_opt("d", "description", "Set image description");
178  Optionpk<bool> align_opt("align", "align", "Align output bounding box to input image",false);
179  Optionpk<double> scale_opt("scale", "scale", "output=scale*input+offset");
180  Optionpk<double> offset_opt("offset", "offset", "output=scale*input+offset");
181  Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
182 
183  extent_opt.setHide(1);
184  cut_opt.setHide(1);
185  eoption_opt.setHide(1);
186  mask_opt.setHide(1);
187  msknodata_opt.setHide(1);
188  mskband_opt.setHide(1);
189  option_opt.setHide(1);
190  file_opt.setHide(1);
191  weight_opt.setHide(1);
192  class_opt.setHide(1);
193  colorTable_opt.setHide(1);
194  description_opt.setHide(1);
195  scale_opt.setHide(1);
196  offset_opt.setHide(1);
197 
198  bool doProcess;//stop process when program was invoked with help option (-h --help)
199  try{
200  doProcess=input_opt.retrieveOption(argc,argv);
201  output_opt.retrieveOption(argc,argv);
202  band_opt.retrieveOption(argc,argv);
203  dx_opt.retrieveOption(argc,argv);
204  dy_opt.retrieveOption(argc,argv);
205  extent_opt.retrieveOption(argc,argv);
206  cut_opt.retrieveOption(argc,argv);
207  eoption_opt.retrieveOption(argc,argv);
208  mask_opt.retrieveOption(argc,argv);
209  msknodata_opt.retrieveOption(argc,argv);
210  mskband_opt.retrieveOption(argc,argv);
211  ulx_opt.retrieveOption(argc,argv);
212  uly_opt.retrieveOption(argc,argv);
213  lrx_opt.retrieveOption(argc,argv);
214  lry_opt.retrieveOption(argc,argv);
215  crule_opt.retrieveOption(argc,argv);
216  ruleBand_opt.retrieveOption(argc,argv);
217  srcnodata_opt.retrieveOption(argc,argv);
218  bndnodata_opt.retrieveOption(argc,argv);
219  minValue_opt.retrieveOption(argc,argv);
220  maxValue_opt.retrieveOption(argc,argv);
221  dstnodata_opt.retrieveOption(argc,argv);
222  resample_opt.retrieveOption(argc,argv);
223  otype_opt.retrieveOption(argc,argv);
224  oformat_opt.retrieveOption(argc,argv);
225  option_opt.retrieveOption(argc,argv);
226  projection_opt.retrieveOption(argc,argv);
227  file_opt.retrieveOption(argc,argv);
228  weight_opt.retrieveOption(argc,argv);
229  class_opt.retrieveOption(argc,argv);
230  colorTable_opt.retrieveOption(argc,argv);
231  description_opt.retrieveOption(argc,argv);
232  align_opt.retrieveOption(argc,argv);
233  scale_opt.retrieveOption(argc,argv);
234  offset_opt.retrieveOption(argc,argv);
235  verbose_opt.retrieveOption(argc,argv);
236  }
237  catch(string predefinedString){
238  std::cout << predefinedString << std::endl;
239  exit(0);
240  }
241  if(!doProcess){
242  cout << endl;
243  cout << "Usage: pkcomposite -i input [-i input]* -o output" << endl;
244  cout << endl;
245  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
246  exit(0);//help was invoked, stop processing
247  }
248 
249  std::map<std::string, crule::CRULE_TYPE> cruleMap;
250  // //initialize cruleMap
251  // enum CRULE_TYPE {overwrite=0, maxndvi=1, maxband=2, minband=3, validband=4, mean=5, mode=6, median=7,sum=8};
252 
253  cruleMap["overwrite"]=crule::overwrite;
254  cruleMap["maxndvi"]=crule::maxndvi;
255  cruleMap["maxband"]=crule::maxband;
256  cruleMap["minband"]=crule::minband;
257  cruleMap["validband"]=crule::validband;
258  cruleMap["mean"]=crule::mean;
259  cruleMap["mode"]=crule::mode;
260  cruleMap["median"]=crule::median;
261  cruleMap["sum"]=crule::sum;
262  cruleMap["maxallbands"]=crule::maxallbands;
263  cruleMap["minallbands"]=crule::minallbands;
264  cruleMap["stdev"]=crule::stdev;
265 
266  if(srcnodata_opt.size()){
267  while(srcnodata_opt.size()<bndnodata_opt.size())
268  srcnodata_opt.push_back(srcnodata_opt[0]);
269  }
270  while(bndnodata_opt.size()<srcnodata_opt.size())
271  bndnodata_opt.push_back(bndnodata_opt[0]);
272  if(minValue_opt.size()){
273  while(minValue_opt.size()<bndnodata_opt.size())
274  minValue_opt.push_back(minValue_opt[0]);
275  while(bndnodata_opt.size()<minValue_opt.size())
276  bndnodata_opt.push_back(bndnodata_opt[0]);
277  }
278  if(maxValue_opt.size()){
279  while(maxValue_opt.size()<bndnodata_opt.size())
280  maxValue_opt.push_back(maxValue_opt[0]);
281  while(bndnodata_opt.size()<maxValue_opt.size())
282  bndnodata_opt.push_back(bndnodata_opt[0]);
283  }
284 
285  RESAMPLE theResample;
286  if(resample_opt[0]=="near"){
287  theResample=NEAR;
288  if(verbose_opt[0])
289  cout << "resampling: nearest neighbor" << endl;
290  }
291  else if(resample_opt[0]=="bilinear"){
292  theResample=BILINEAR;
293  if(verbose_opt[0])
294  cout << "resampling: bilinear interpolation" << endl;
295  }
296  else{
297  std::cout << "Error: resampling method " << resample_opt[0] << " not supported" << std::endl;
298  exit(1);
299  }
300 
301  if(input_opt.empty()){
302  std::cerr << "No input file provided (use option -i). Use --help for help information" << std::endl;
303  exit(0);
304  }
305  int nband=0;
306  int nwriteBand=0;
307  int writeBand=0;
308  vector<short> bands;
309 
310  //get bounding box
311  double maxLRX=lrx_opt[0];
312  double maxULY=uly_opt[0];
313  double minULX=ulx_opt[0];
314  double minLRY=lry_opt[0];
315  double magic_x=1,magic_y=1;//magic pixel for GDAL map info
316 
317  GDALDataType theType=GDT_Unknown;
318  if(verbose_opt[0])
319  cout << "possible output data types: ";
320  for(int iType = 0; iType < GDT_TypeCount; ++iType){
321  if(verbose_opt[0])
322  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
323  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
324  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
325  otype_opt[0].c_str()))
326  theType=(GDALDataType) iType;
327  }
328  if(verbose_opt[0]){
329  cout << endl;
330  if(theType==GDT_Unknown)
331  cout << "Unknown output pixel type: " << otype_opt[0] << endl;
332  else
333  cout << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
334  }
335 
336  double dx=0;
337  double dy=0;
338  //get bounding box from extentReader if defined
339  ImgReaderOgr extentReader;
340  if(extent_opt.size()){
341  double e_ulx;
342  double e_uly;
343  double e_lrx;
344  double e_lry;
345  for(int iextent=0;iextent<extent_opt.size();++iextent){
346  extentReader.open(extent_opt[iextent]);
347  if(!(extentReader.getExtent(e_ulx,e_uly,e_lrx,e_lry))){
348  cerr << "Error: could not get extent from " << extent_opt[0] << endl;
349  exit(1);
350  }
351  if(!iextent){
352  ulx_opt[0]=e_ulx;
353  uly_opt[0]=e_uly;
354  lrx_opt[0]=e_lrx;
355  lry_opt[0]=e_lry;
356  }
357  else{
358  if(e_ulx<ulx_opt[0])
359  ulx_opt[0]=e_ulx;
360  if(e_uly>uly_opt[0])
361  uly_opt[0]=e_uly;
362  if(e_lrx>lrx_opt[0])
363  lrx_opt[0]=e_lrx;
364  if(e_lry<lry_opt[0])
365  lry_opt[0]=e_lry;
366  }
367  extentReader.close();
368  }
369  if(maxLRX>minULX&&minULX>ulx_opt[0])
370  ulx_opt[0]=minULX;
371  if(maxLRX>minULX&&maxLRX<lrx_opt[0])
372  lrx_opt[0]=maxLRX;
373  if(maxULY>minLRY&&maxULY<uly_opt[0])
374  uly_opt[0]=maxULY;
375  if(minLRY<maxULY&&minLRY>lry_opt[0])
376  lry_opt[0]=minLRY;
377  if(cut_opt.size()||eoption_opt.size())
378  extentReader.open(extent_opt[0]);
379  }
380 
381  if(verbose_opt[0])
382  cout << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
383 
384  vector<ImgReaderGdal> imgReader(input_opt.size());
385  // vector<ImgReaderGdal> imgReader(input_opt.size());
386  string theProjection="";
387  GDALColorTable* theColorTable=NULL;
388  string imageType;
389  bool init=false;
390  for(int ifile=0;ifile<input_opt.size();++ifile){
391  try{
392  imgReader[ifile].open(input_opt[ifile],GA_ReadOnly);
393  // imgReader[ifile].open(input_opt[ifile],GA_ReadOnly);
394  for(int iband=0;iband<scale_opt.size();++iband)
395  imgReader[ifile].setScale(scale_opt[iband],iband);
396  for(int iband=0;iband<offset_opt.size();++iband)
397  imgReader[ifile].setOffset(offset_opt[iband],iband);
398  }
399  catch(string errorstring){
400  cerr << errorstring << " " << input_opt[ifile] << endl;
401  }
402  //todo: must be in init part only?
403  if(colorTable_opt.empty())
404  if(imgReader[ifile].getColorTable())
405  theColorTable=(imgReader[ifile].getColorTable()->Clone());
406  if(projection_opt.empty())
407  theProjection=imgReader[ifile].getProjection();
408  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
409  string theInterleave="INTERLEAVE=";
410  theInterleave+=imgReader[ifile].getInterleave();
411  option_opt.push_back(theInterleave);
412  }
413 
414  if((ulx_opt[0]||uly_opt[0]||lrx_opt[0]||lry_opt[0])&&(!imgReader[ifile].covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
415  if(verbose_opt[0])
416  cout << input_opt[ifile] << " not within bounding box, skipping..." << endl;
417  continue;
418  }
419  double theULX, theULY, theLRX, theLRY;
420  imgReader[ifile].getBoundingBox(theULX,theULY,theLRX,theLRY);
421  if(theLRY>theULY){
422  cerr << "Error: " << input_opt[ifile] << " is not georeferenced, only referenced images are supported for pkcomposite " << endl;
423  exit(1);
424  }
425  if(verbose_opt[0])
426  cout << "Bounding Box (ULX ULY LRX LRY): " << fixed << setprecision(6) << theULX << " " << theULY << " " << theLRX << " " << theLRY << endl;
427  if(!init){
428  if(verbose_opt[0]){
429  switch(cruleMap[crule_opt[0]]){
430  default:
431  case(crule::overwrite):
432  cout << "Composite rule: overwrite" << endl;
433  break;
434  case(crule::maxndvi):
435  cout << "Composite rule: max ndvi" << endl;
436  break;
437  case(crule::maxband):
438  cout << "Composite rule: max band" << endl;
439  break;
440  case(crule::minband):
441  cout << "Composite rule: min band" << endl;
442  break;
443  case(crule::validband):
444  cout << "Composite rule: valid band" << endl;
445  break;
446  case(crule::mean):
447  cout << "Composite rule: mean value" << endl;
448  break;
449  case(crule::mode):
450  cout << "Composite rule: max voting (only for byte images)" << endl;
451  break;
452  case(crule::median):
453  cout << "Composite rule: median" << endl;
454  break;
455  case(crule::stdev):
456  cout << "Composite rule: stdev" << endl;
457  break;
458  case(crule::sum):
459  cout << "Composite rule: sum" << endl;
460  break;
461  case(crule::minallbands):
462  cout << "Composite rule: minallbands" << endl;
463  break;
464  case(crule::maxallbands):
465  cout << "Composite rule: maxallbands" << endl;
466  break;
467  }
468  }
469  if(band_opt.size()){
470  nband=band_opt.size();
471  bands.resize(band_opt.size());
472  for(int iband=0;iband<band_opt.size();++iband){
473  bands[iband]=band_opt[iband];
474  assert(bands[iband]<imgReader[ifile].nrOfBand());
475  }
476  }
477  else{
478  nband=imgReader[ifile].nrOfBand();
479  bands.resize(nband);
480  for(int iband=0;iband<nband;++iband)
481  bands[iband]=iband;
482  }
483  for(int iband=0;iband<bndnodata_opt.size();++iband){
484  assert(bndnodata_opt[iband]>=0&&bndnodata_opt[iband]<nband);
485  }
486  //if output type not set, get type from input image
487  if(theType==GDT_Unknown){
488  theType=imgReader[ifile].getDataType();
489  if(verbose_opt[0])
490  cout << "Using data type from input image: " << GDALGetDataTypeName(theType) << endl;
491  }
492 
493  if(oformat_opt.size())//default
494  imageType=oformat_opt[0];
495  else
496  imageType=imgReader[ifile].getImageType();
497 
498  if(verbose_opt[0]){
499  cout << "type of data for " << input_opt[ifile] << ": " << theType << endl;
500  cout << "nband: " << nband << endl;
501  }
502 
503  maxLRX=theLRX;
504  maxULY=theULY;
505  minULX=theULX;
506  minLRY=theLRY;
507  if(dx_opt.size())
508  dx=dx_opt[0];
509  else
510  dx=imgReader[ifile].getDeltaX();
511  if(dy_opt.size())
512  dy=dy_opt[0];
513  else
514  dy=imgReader[ifile].getDeltaY();
515  init=true;
516  }
517  else{
518  maxLRX=(theLRX>maxLRX)?theLRX:maxLRX;
519  maxULY=(theULY>maxULY)?theULY:maxULY;
520  minULX=(theULX<minULX)?theULX:minULX;
521  minLRY=(theLRY<minLRY)?theLRY:minLRY;
522  }
523  // imgReader.close();
524  }
525  if(verbose_opt[0])
526  cout << "bounding box input images (ULX ULY LRX LRY): " << fixed << setprecision(6) << minULX << " " << maxULY << " " << maxLRX << " " << minLRY << endl;
527  if(ulx_opt[0]||uly_opt[0]||lrx_opt[0]||lry_opt[0]){
528  maxLRX=lrx_opt[0];
529  maxULY=uly_opt[0];
530  minULX=ulx_opt[0];
531  minLRY=lry_opt[0];
532  }
533 
534  bool forceEUgrid=false;
535  if(projection_opt.size())
536  forceEUgrid=(!(projection_opt[0].compare("EPSG:3035"))||!(projection_opt[0].compare("EPSG:3035"))||projection_opt[0].find("ETRS-LAEA")!=string::npos);
537  if(forceEUgrid){
538  //force to LAEA grid
539  minULX=floor(minULX);
540  minULX-=static_cast<int>(minULX)%(static_cast<int>(dx));
541  maxULY=ceil(maxULY);
542  if(static_cast<int>(maxULY)%static_cast<int>(dy))
543  maxULY+=dy;
544  maxULY-=static_cast<int>(maxULY)%(static_cast<int>(dy));
545  maxLRX=ceil(maxLRX);
546  if(static_cast<int>(maxLRX)%static_cast<int>(dx))
547  maxLRX+=dx;
548  maxLRX-=static_cast<int>(maxLRX)%(static_cast<int>(dx));
549  minLRY=floor(minLRY);
550  minLRY-=static_cast<int>(minLRY)%(static_cast<int>(dy));
551  }
552  else if(align_opt[0]){
553  if(minULX>imgReader[0].getUlx())
554  minULX-=fmod(minULX-imgReader[0].getUlx(),dx);
555  else if(minULX<imgReader[0].getUlx())
556  minULX+=fmod(imgReader[0].getUlx()-minULX,dx)-dx;
557  if(maxLRX<imgReader[0].getLrx())
558  maxLRX+=fmod(imgReader[0].getLrx()-maxLRX,dx);
559  else if(maxLRX>imgReader[0].getLrx())
560  maxLRX-=fmod(maxLRX-imgReader[0].getLrx(),dx)+dx;
561  if(minLRY>imgReader[0].getLry())
562  minLRY-=fmod(minLRY-imgReader[0].getLry(),dy);
563  else if(minLRY<imgReader[0].getLry())
564  minLRY+=fmod(imgReader[0].getLry()-minLRY,dy)-dy;
565  if(maxULY<imgReader[0].getUly())
566  maxULY+=fmod(imgReader[0].getUly()-maxULY,dy);
567  else if(maxULY>imgReader[0].getUly())
568  maxULY-=fmod(maxULY-imgReader[0].getUly(),dy)+dy;
569  }
570 
571  if(verbose_opt[0])
572  cout << "bounding box composite image (ULX ULY LRX LRY): " << fixed << setprecision(6) << minULX << " " << maxULY << " " << maxLRX << " " << minLRY << endl;
573  //initialize image
574  if(verbose_opt[0])
575  cout << "initializing composite image..." << endl;
576 // double dcol=(maxLRX-minULX+dx-1)/dx;
577 // double drow=(maxULY-minLRY+dy-1)/dy;
578 // int ncol=static_cast<int>(dcol);
579 // int nrow=static_cast<int>(drow);
580 
581  int ncol=ceil((maxLRX-minULX)/dx);
582  int nrow=ceil((maxULY-minLRY)/dy);
583 
584  if(verbose_opt[0])
585  cout << "composite image dim (nrow x ncol): " << nrow << " x " << ncol << endl;
586  ImgWriterGdal imgWriter;
587  while(weight_opt.size()<input_opt.size())
588  weight_opt.push_back(weight_opt[0]);
589  if(verbose_opt[0]){
590  std::cout << weight_opt << std::endl;
591  }
592  if(cruleMap[crule_opt[0]]==crule::mode){
593  nwriteBand=(file_opt[0])? class_opt.size()+1:class_opt.size();
594  }
595  else
596  nwriteBand=(file_opt[0])? bands.size()+1:bands.size();
597  if(output_opt.empty()){
598  std::cerr << "No output file provided (use option -o). Use --help for help information" << std::endl;
599  exit(0);
600  }
601  if(verbose_opt[0])
602  cout << "open output image " << output_opt[0] << " with " << nwriteBand << " bands" << endl << flush;
603 
604  try{
605  imgWriter.open(output_opt[0],ncol,nrow,nwriteBand,theType,imageType,option_opt);
606  imgWriter.setNoData(dstnodata_opt);
607  }
608  catch(string error){
609  cout << error << endl;
610  }
611  if(description_opt.size())
612  imgWriter.setImageDescription(description_opt[0]);
613  double gt[6];
614  gt[0]=minULX;
615  gt[1]=dx;
616  gt[2]=0;
617  gt[3]=maxULY;
618  gt[4]=0;
619  gt[5]=-dy;
620  imgWriter.setGeoTransform(gt);
621  // imgWriter.setGeoTransform(minULX,maxULY,dx,dy,0,0);
622  if(projection_opt.size()){
623  if(verbose_opt[0])
624  cout << "projection: " << projection_opt[0] << endl;
625  imgWriter.setProjectionProj4(projection_opt[0]);
626  }
627  else if(theProjection!=""){
628  if(verbose_opt[0])
629  cout << "projection: " << theProjection << endl;
630  imgWriter.setProjection(theProjection);
631  }
632  if(imgWriter.getDataType()==GDT_Byte){
633  if(colorTable_opt.size()){
634  if(colorTable_opt[0]!="none")
635  imgWriter.setColorTable(colorTable_opt[0]);
636  }
637  else if(theColorTable)
638  imgWriter.setColorTable(theColorTable);
639  }
640 
641  ImgWriterGdal maskWriter;
642  if(extent_opt.size()&&(cut_opt[0]||eoption_opt.size())){
643  try{
644  maskWriter.open("/vsimem/mask.tif",ncol,nrow,1,GDT_Float32,"GTiff",option_opt);
645  double gt[6];
646  gt[0]=minULX;
647  gt[1]=dx;
648  gt[2]=0;
649  gt[3]=maxULY;
650  gt[4]=0;
651  gt[5]=-dy;
652  maskWriter.setGeoTransform(gt);
653  if(projection_opt.size())
654  maskWriter.setProjectionProj4(projection_opt[0]);
655  else if(theProjection!=""){
656  if(verbose_opt[0])
657  cout << "projection: " << theProjection << endl;
658  maskWriter.setProjection(theProjection);
659  }
660  vector<double> burnValues(1,1);//burn value is 1 (single band)
661  maskWriter.rasterizeOgr(extentReader,burnValues,eoption_opt);
662  maskWriter.close();
663  }
664  catch(string error){
665  cerr << error << std::endl;
666  exit(2);
667  }
668  catch(...){
669  cerr << "error caught" << std::endl;
670  exit(1);
671  }
672  //todo: support multiple masks
673  mask_opt.clear();
674  mask_opt.push_back("/vsimem/mask.tif");
675  }
676  ImgReaderGdal maskReader;
677  if(mask_opt.size()){
678  try{
679  if(verbose_opt[0]>=1)
680  std::cout << "opening mask image file " << mask_opt[0] << std::endl;
681  maskReader.open(mask_opt[0],GA_ReadOnly);
682  if(mskband_opt[0]>=maskReader.nrOfBand()){
683  string errorString="Error: illegal mask band";
684  throw(errorString);
685  }
686  }
687  catch(string error){
688  cerr << error << std::endl;
689  exit(2);
690  }
691  catch(...){
692  cerr << "error caught" << std::endl;
693  exit(1);
694  }
695  }
696 
697  //create composite image
698  if(verbose_opt[0])
699  cout << "creating composite image" << endl;
700  Vector2d<double> writeBuffer(nband,imgWriter.nrOfCol());
701  vector<short> fileBuffer(ncol);//holds the number of used files
702  Vector2d<short> maxBuffer;//buffer used for maximum voting
703  // Vector2d<double> readBuffer(nband);
704  vector<Vector2d<unsigned short> > readBuffer(input_opt.size());
705  for(int ifile=0;ifile<input_opt.size();++ifile)
706  readBuffer[ifile].resize(imgReader[ifile].nrOfBand());
708  if(cruleMap[crule_opt[0]]==crule::maxndvi)//ndvi
709  assert(ruleBand_opt.size()==2);
710  if(cruleMap[crule_opt[0]]==crule::mode){//max voting
711  maxBuffer.resize(imgWriter.nrOfCol(),256);//use only byte images for max voting
712  for(int iclass=0;iclass<class_opt.size();++iclass)
713  assert(class_opt[iclass]<maxBuffer.size());
714  }
715  int jb=0;
716  double readRow=0;
717  double readCol=0;
718  double lowerCol=0;
719  double upperCol=0;
720  const char* pszMessage;
721  void* pProgressArg=NULL;
722  GDALProgressFunc pfnProgress=GDALTermProgress;
723  double progress=0;
724  pfnProgress(progress,pszMessage,pProgressArg);
725  for(int irow=0;irow<imgWriter.nrOfRow();++irow){
726  vector<float> lineMask;
727  Vector2d< vector<double> > storeBuffer;
728  vector<bool> writeValid(ncol);
729 
730  //convert irow to geo
731  double x=0;
732  double y=0;
733  imgWriter.image2geo(0,irow,x,y);
734 
735 
736  if(cruleMap[crule_opt[0]]==crule::mean ||
737  cruleMap[crule_opt[0]]==crule::median ||
738  cruleMap[crule_opt[0]]==crule::sum ||
739  cruleMap[crule_opt[0]]==crule::minallbands ||
740  cruleMap[crule_opt[0]]==crule::maxallbands ||
741  cruleMap[crule_opt[0]]==crule::stdev)
742  storeBuffer.resize(nband,ncol);
743  for(int icol=0;icol<imgWriter.nrOfCol();++icol){
744  writeValid[icol]=false;
745  fileBuffer[icol]=0;
746  if(cruleMap[crule_opt[0]]==crule::mode){//max voting
747  for(int iclass=0;iclass<256;++iclass)
748  maxBuffer[icol][iclass]=0;
749  }
750  else{
751  for(int iband=0;iband<nband;++iband)
752  writeBuffer[iband][icol]=dstnodata_opt[0];
753  }
754  }
755 
756  double oldRowMask=-1;//keep track of row mask to optimize number of line readings
757 
758  for(int ifile=0;ifile<input_opt.size();++ifile){
759  //imgReader already open...
760  // try{
761  // imgReader.open(input_opt[ifile]);
762  // }
763  // catch(string error){
764  // cout << error << endl;
765  // }
766  // assert(imgReader.getDataType()==theType);
767  assert(imgReader[ifile].nrOfBand()>=nband);
768  if(!imgReader[ifile].covers(minULX,maxULY,maxLRX,minLRY)){
769  // imgReader.close();
770  continue;
771  }
772  double uli,ulj,lri,lrj;
773  imgReader[ifile].geo2image(minULX+(magic_x-1.0)*imgReader[ifile].getDeltaX(),maxULY-(magic_y-1.0)*imgReader[ifile].getDeltaY(),uli,ulj);
774  imgReader[ifile].geo2image(maxLRX+(magic_x-2.0)*imgReader[ifile].getDeltaX(),minLRY-(magic_y-2.0)*imgReader[ifile].getDeltaY(),lri,lrj);
775  uli=floor(uli);
776  ulj=floor(ulj);
777  lri=floor(lri);
778  lrj=floor(lrj);
779 
780  double startCol=uli;
781  double endCol=lri;
782  if(uli<0)
783  startCol=0;
784  else if(uli>=imgReader[ifile].nrOfCol())
785  startCol=imgReader[ifile].nrOfCol()-1;
786  if(lri<0)
787  endCol=0;
788  else if(lri>=imgReader[ifile].nrOfCol())
789  endCol=imgReader[ifile].nrOfCol()-1;
790  int readncol=endCol-startCol+1;
791 
792  //lookup corresponding row for irow in this file
793  imgReader[ifile].geo2image(x,y,readCol,readRow);
794  if(readRow<0||readRow>=imgReader[ifile].nrOfRow()){
795  // imgReader.close();
796  continue;
797  }
798  // for(int iband=0;iband<imgReader.nrOfBand();++iband){
799  for(int iband=0;iband<nband;++iband){
800  int readBand=(band_opt.size()>iband)? band_opt[iband] : iband;
801  // readBuffer[iband].resize(readncol);
802  try{
803  imgReader[ifile].readData(readBuffer[ifile][iband],startCol,endCol,readRow,readBand,theResample);
804  //test
805  // imgReader[ifile].readData(readBuffer[ifile][iband],static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(readRow),static_cast<int>(readBand));
806  // if(readRow==0&&iband==0){
807  // for(int icol=0;icol<10;++icol)
808  // cout << readBuffer[0][0][icol] << " ";
809  // cout << endl;
810  // }
811  }
812  catch(string error){
813  cerr << "error reading image " << input_opt[ifile] << ": " << endl;
814  throw;
815  }
816  }
817  for(int ib=0;ib<ncol;++ib){
818  imgWriter.image2geo(ib,irow,x,y);
819  //check mask first
820  bool valid=true;
821  if(mask_opt.size()){
822  //read mask
823  double colMask=0;
824  double rowMask=0;
825 
826  maskReader.geo2image(x,y,colMask,rowMask);
827  colMask=static_cast<int>(colMask);
828  rowMask=static_cast<int>(rowMask);
829  if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
830  if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
831 
832  assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
833  try{
834  maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
835  }
836  catch(string errorstring){
837  cerr << errorstring << endl;
838  exit(1);
839  }
840  catch(...){
841  cerr << "error caught" << std::endl;
842  exit(3);
843  }
844  oldRowMask=rowMask;
845  }
846  if(lineMask[colMask]==msknodata_opt[0])
847  valid=false;
848  }
849  }
850 
851  if(!valid)
852  continue;
853 
854  //lookup corresponding row for irow in this file
855  imgReader[ifile].geo2image(x,y,readCol,readRow);
856  if(readCol<0||readCol>=imgReader[ifile].nrOfCol())
857  continue;
858  double val_current=0;
859  double val_new=0;
860  bool readValid=true;
861  switch(theResample){
862  case(BILINEAR):
863  lowerCol=readCol-0.5;
864  lowerCol=static_cast<int>(lowerCol);
865  upperCol=readCol+0.5;
866  upperCol=static_cast<int>(upperCol);
867  if(lowerCol<0)
868  lowerCol=0;
869  if(upperCol>=imgReader[ifile].nrOfCol())
870  upperCol=imgReader[ifile].nrOfCol()-1;
871  for(int vband=0;vband<bndnodata_opt.size();++vband){
872  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][bndnodata_opt[vband]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][bndnodata_opt[vband]][lowerCol-startCol];
873  if(minValue_opt.size()>vband){
874  if(val_new<=minValue_opt[vband]){
875  readValid=false;
876  break;
877  }
878  }
879  if(maxValue_opt.size()>vband){
880  if(val_new>=maxValue_opt[vband]){
881  readValid=false;
882  break;
883  }
884  }
885  if(srcnodata_opt.size()>vband){
886  if(val_new==srcnodata_opt[vband]){
887  readValid=false;
888  break;
889  }
890  }
891  }
892  break;
893  default:
894  readCol=static_cast<int>(readCol);
895  for(int vband=0;vband<bndnodata_opt.size();++vband){
896  val_new=readBuffer[ifile][bndnodata_opt[vband]][readCol-startCol];
897  if(minValue_opt.size()>vband){
898  if(val_new<=minValue_opt[vband]){
899  readValid=false;
900  break;
901  }
902  }
903  if(maxValue_opt.size()>vband){
904  if(val_new>=maxValue_opt[vband]){
905  readValid=false;
906  break;
907  }
908  }
909  if(srcnodata_opt.size()>vband){
910  if(val_new==srcnodata_opt[vband]){
911  readValid=false;
912  break;
913  }
914  }
915  }
916  break;
917  }
918  if(readValid){
919  if(file_opt[0]==1)
920  ++fileBuffer[ib];
921  if(writeValid[ib]){
922  int iband=0;
923  switch(cruleMap[crule_opt[0]]){
924  case(crule::maxndvi):{//max ndvi
925  double red_current=writeBuffer[ruleBand_opt[0]][ib];
926  double nir_current=writeBuffer[ruleBand_opt[1]][ib];
927  double ndvi_current=0;
928  if(red_current+nir_current>0&&red_current>=0&&nir_current>=0)
929  ndvi_current=(nir_current-red_current)/(nir_current+red_current);
930  double ndvi_new=0;
931  double red_new=0;
932  double nir_new=0;
933  switch(theResample){
934  case(BILINEAR):
935  lowerCol=readCol-0.5;
936  lowerCol=static_cast<int>(lowerCol);
937  upperCol=readCol+0.5;
938  upperCol=static_cast<int>(upperCol);
939  if(lowerCol<0)
940  lowerCol=0;
941  if(upperCol>=imgReader[ifile].nrOfCol())
942  upperCol=imgReader[ifile].nrOfCol()-1;
943  red_new=(readCol-0.5-lowerCol)*readBuffer[ifile][ruleBand_opt[0]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][ruleBand_opt[0]][lowerCol-startCol];
944  nir_new=(readCol-0.5-lowerCol)*readBuffer[ifile][ruleBand_opt[1]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][ruleBand_opt[1]][lowerCol-startCol];
945  if(red_new+nir_new>0&&red_new>=0&&nir_new>=0)
946  ndvi_new=(nir_new-red_new)/(nir_new+red_new);
947  if(ndvi_new>=ndvi_current){
948  for(iband=0;iband<nband;++iband){
949  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
950  writeBuffer[iband][ib]=val_new;
951  }
952  if(file_opt[0]>1)
953  fileBuffer[ib]=ifile;
954  }
955  break;
956  default:
957  readCol=static_cast<int>(readCol);
958  red_new=readBuffer[ifile][ruleBand_opt[0]][readCol-startCol];
959  nir_new=readBuffer[ifile][ruleBand_opt[1]][readCol-startCol];
960  if(red_new+nir_new>0&&red_new>=0&&nir_new>=0)
961  ndvi_new=(nir_new-red_new)/(nir_new+red_new);
962  if(ndvi_new>=ndvi_current){
963  for(iband=0;iband<nband;++iband){
964  val_new=readBuffer[ifile][iband][readCol-startCol];
965  writeBuffer[iband][ib]=val_new;
966  }
967  if(file_opt[0]>1)
968  fileBuffer[ib]=ifile;
969  }
970  break;
971  }
972  break;
973  }
974  case(crule::maxband):
975  case(crule::minband):
976  case(crule::validband)://max,min,valid band
977  val_current=writeBuffer[ruleBand_opt[0]][ib];
978  switch(theResample){
979  case(BILINEAR):
980  lowerCol=readCol-0.5;
981  lowerCol=static_cast<int>(lowerCol);
982  upperCol=readCol+0.5;
983  upperCol=static_cast<int>(upperCol);
984  if(lowerCol<0)
985  lowerCol=0;
986  if(upperCol>=imgReader[ifile].nrOfCol())
987  upperCol=imgReader[ifile].nrOfCol()-1;
988  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][ruleBand_opt[0]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][ruleBand_opt[0]][lowerCol-startCol];
989  val_new*=weight_opt[ifile];
990  if((cruleMap[crule_opt[0]]==crule::maxband&&val_new>val_current)||(cruleMap[crule_opt[0]]==crule::minband&&val_new<val_current)||(cruleMap[crule_opt[0]]==crule::validband)){//&&val_new>minValue_opt[0]&&val_new<maxValue_opt[0])){
991  for(iband=0;iband<nband;++iband){
992  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
993  val_new*=weight_opt[ifile];
994  writeBuffer[iband][ib]=val_new;
995  }
996  if(file_opt[0]>1)
997  fileBuffer[ib]=ifile;
998  }
999  break;
1000  default:
1001  readCol=static_cast<int>(readCol);
1002  val_new=readBuffer[ifile][ruleBand_opt[0]][readCol-startCol];
1003  val_new*=weight_opt[ifile];
1004  if((cruleMap[crule_opt[0]]==crule::maxband&&val_new>val_current)||(cruleMap[crule_opt[0]]==crule::minband&&val_new<val_current)||(cruleMap[crule_opt[0]]==crule::validband)){//&&val_new>minValue_opt[0]&&val_new<maxValue_opt[0])){
1005  for(iband=0;iband<nband;++iband){
1006  val_new=readBuffer[ifile][iband][readCol-startCol];
1007  val_new*=weight_opt[ifile];
1008  writeBuffer[iband][ib]=val_new;
1009  }
1010  if(file_opt[0]>1)
1011  fileBuffer[ib]=ifile;
1012  }
1013  break;
1014  }
1015  break;
1016  case(crule::mode)://max voting (only for Byte images)
1017  switch(theResample){
1018  case(BILINEAR):
1019  lowerCol=readCol-0.5;
1020  lowerCol=static_cast<int>(lowerCol);
1021  upperCol=readCol+0.5;
1022  upperCol=static_cast<int>(upperCol);
1023  if(lowerCol<0)
1024  lowerCol=0;
1025  if(upperCol>=imgReader[ifile].nrOfCol())
1026  upperCol=imgReader[ifile].nrOfCol()-1;
1027  for(iband=0;iband<nband;++iband){
1028  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
1029  maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
1030  // ++(maxBuffer[ib][val_new]);
1031  }
1032  break;
1033  default:
1034  readCol=static_cast<int>(readCol);
1035  for(iband=0;iband<nband;++iband){
1036  val_new=readBuffer[ifile][iband][readCol-startCol];
1037  maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
1038  }
1039  break;
1040  }
1041  break;
1042  case(crule::mean)://mean value
1043  case(crule::median)://median value
1044  case(crule::sum)://sum value
1045  case(crule::minallbands)://minimum for each and every band
1046  case(crule::maxallbands)://maximum for each and every band
1047  case(crule::stdev)://maximum for each and every band
1048  switch(theResample){
1049  case(BILINEAR):
1050  lowerCol=readCol-0.5;
1051  lowerCol=static_cast<int>(lowerCol);
1052  upperCol=readCol+0.5;
1053  upperCol=static_cast<int>(upperCol);
1054  if(lowerCol<0)
1055  lowerCol=0;
1056  if(upperCol>=imgReader[ifile].nrOfCol())
1057  upperCol=imgReader[ifile].nrOfCol()-1;
1058  for(iband=0;iband<nband;++iband){
1059  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
1060  val_new*=weight_opt[ifile];
1061  storeBuffer[iband][ib].push_back(val_new);
1062  }
1063  break;
1064  default:
1065  readCol=static_cast<int>(readCol);
1066  for(iband=0;iband<nband;++iband){
1067  val_new=readBuffer[ifile][iband][readCol-startCol];
1068  val_new*=weight_opt[ifile];
1069  storeBuffer[iband][ib].push_back(val_new);
1070  assert(ifile>0);
1071  // assert(weight_opt[ifile]>=0);
1072  // assert(storeBuffer[iband][ib].back()>=0);
1073  }
1074  break;
1075  }
1076  if(file_opt[0]>1)
1077  fileBuffer[ib]=ifile;
1078  break;
1079  case(crule::overwrite):
1080  default:
1081  switch(theResample){
1082  case(BILINEAR):
1083  lowerCol=readCol-0.5;
1084  lowerCol=static_cast<int>(lowerCol);
1085  upperCol=readCol+0.5;
1086  upperCol=static_cast<int>(upperCol);
1087  if(lowerCol<0)
1088  lowerCol=0;
1089  if(upperCol>=imgReader[ifile].nrOfCol())
1090  upperCol=imgReader[ifile].nrOfCol()-1;
1091  for(iband=0;iband<nband;++iband){
1092  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
1093  val_new*=weight_opt[ifile];
1094  writeBuffer[iband][ib]=val_new;
1095  }
1096  break;
1097  default:
1098  readCol=static_cast<int>(readCol);
1099  for(iband=0;iband<nband;++iband){
1100  val_new=readBuffer[ifile][iband][readCol-startCol];
1101  val_new*=weight_opt[ifile];
1102  writeBuffer[iband][ib]=val_new;
1103  }
1104  break;
1105  }
1106  if(file_opt[0]>1)
1107  fileBuffer[ib]=ifile;
1108  break;
1109  }
1110  }
1111  else{
1112  writeValid[ib]=true;//readValid was true
1113  int iband=0;
1114  switch(cruleMap[crule_opt[0]]){
1115  case(crule::mean):
1116  case(crule::median):
1117  case(crule::sum):
1118  case(crule::minallbands):
1119  case(crule::maxallbands):
1120  case(crule::stdev):
1121  switch(theResample){
1122  case(BILINEAR):
1123  lowerCol=readCol-0.5;
1124  lowerCol=static_cast<int>(lowerCol);
1125  upperCol=readCol+0.5;
1126  upperCol=static_cast<int>(upperCol);
1127  if(lowerCol<0)
1128  lowerCol=0;
1129  if(upperCol>=imgReader[ifile].nrOfCol())
1130  upperCol=imgReader[ifile].nrOfCol()-1;
1131  for(iband=0;iband<nband;++iband){
1132  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
1133  val_new*=weight_opt[ifile];
1134  storeBuffer[iband][ib].push_back(val_new);
1135  }
1136  break;
1137  default:
1138  readCol=static_cast<int>(readCol);
1139  for(iband=0;iband<nband;++iband){
1140  val_new=readBuffer[ifile][iband][readCol-startCol];
1141  val_new*=weight_opt[ifile];
1142  storeBuffer[iband][ib].push_back(val_new);
1143  }
1144  break;
1145  }
1146  if(file_opt[0]>1)
1147  fileBuffer[ib]=ifile;
1148  break;
1149  case(crule::mode):
1150  switch(theResample){
1151  case(BILINEAR):
1152  lowerCol=readCol-0.5;
1153  lowerCol=static_cast<int>(lowerCol);
1154  upperCol=readCol+0.5;
1155  upperCol=static_cast<int>(upperCol);
1156  if(lowerCol<0)
1157  lowerCol=0;
1158  if(upperCol>=imgReader[ifile].nrOfCol())
1159  upperCol=imgReader[ifile].nrOfCol()-1;
1160  for(iband=0;iband<nband;++iband){
1161  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
1162  maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
1163  // ++(maxBuffer[ib][val_new]);
1164  }
1165  break;
1166  default:
1167  readCol=static_cast<int>(readCol);
1168  for(iband=0;iband<nband;++iband){
1169  val_new=readBuffer[ifile][iband][readCol-startCol];
1170  maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
1171  }
1172  // ++(maxBuffer[ib][val_new]);
1173  break;
1174  }
1175  break;
1176  default:
1177  switch(theResample){
1178  case(BILINEAR):
1179  lowerCol=readCol-0.5;
1180  lowerCol=static_cast<int>(lowerCol);
1181  upperCol=readCol+0.5;
1182  upperCol=static_cast<int>(upperCol);
1183  if(lowerCol<0)
1184  lowerCol=0;
1185  if(upperCol>=imgReader[ifile].nrOfCol())
1186  upperCol=imgReader[ifile].nrOfCol()-1;
1187  for(iband=0;iband<nband;++iband){
1188  val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
1189  val_new*=weight_opt[ifile];
1190  writeBuffer[iband][ib]=val_new;
1191  }
1192  break;
1193  default:
1194  readCol=static_cast<int>(readCol);
1195  for(iband=0;iband<nband;++iband){
1196  val_new=readBuffer[ifile][iband][readCol-startCol];
1197  val_new*=weight_opt[ifile];
1198  writeBuffer[iband][ib]=val_new;
1199  }
1200  break;
1201  }
1202  if(file_opt[0]>1)
1203  fileBuffer[ib]=ifile;
1204  break;
1205  }
1206  }
1207  }
1208  }
1209  // imgReader.close();
1210  }
1211  if(cruleMap[crule_opt[0]]==crule::mode){
1212  vector<short> classBuffer(imgWriter.nrOfCol());
1213  if(class_opt.size()>1){
1214  for(int iclass=0;iclass<class_opt.size();++iclass){
1215  for(int icol=0;icol<imgWriter.nrOfCol();++icol)
1216  classBuffer[icol]=maxBuffer[icol][class_opt[iclass]];
1217  try{
1218  imgWriter.writeData(classBuffer,irow,iclass);
1219  }
1220  catch(string error){
1221  cerr << "error writing image file " << output_opt[0] << ": " << error << endl;
1222  throw;
1223  }
1224  }
1225  }
1226  else{
1227  for(int icol=0;icol<imgWriter.nrOfCol();++icol){
1228  vector<short>::iterator maxit=maxBuffer[icol].begin();
1229  maxit=stat.mymax(maxBuffer[icol],maxBuffer[icol].begin(),maxBuffer[icol].end());
1230  writeBuffer[0][icol]=distance(maxBuffer[icol].begin(),maxit);
1231  if(file_opt[0]>1)
1232  fileBuffer[icol]=*(maxit);
1233  }
1234  try{
1235  imgWriter.writeData(writeBuffer[0],irow,0);
1236  if(file_opt[0])
1237  imgWriter.writeData(fileBuffer,irow,1);
1238  }
1239  catch(string error){
1240  cerr << "error writing image file " << output_opt[0] << ": " << error << endl;
1241  throw;
1242  }
1243  }
1244  }
1245  else{
1246  for(int iband=0;iband<bands.size();++iband){
1247  // assert(writeBuffer[bands[iband]].size()==imgWriter.nrOfCol());
1248  assert(writeBuffer[iband].size()==imgWriter.nrOfCol());
1249  for(int icol=0;icol<imgWriter.nrOfCol();++icol){
1250  try{
1251  switch(cruleMap[crule_opt[0]]){
1252  case(crule::mean):
1253  // writeBuffer[iband][icol]=stat.mean(storeBuffer[bands[iband]][icol]);
1254  writeBuffer[iband][icol]=stat.mean(storeBuffer[iband][icol]);
1255  break;
1256  case(crule::median):
1257  // writeBuffer[iband][icol]=stat.median(storeBuffer[bands[iband]][icol]);
1258  writeBuffer[iband][icol]=stat.median(storeBuffer[iband][icol]);
1259  break;
1260  case(crule::sum):
1261  // writeBuffer[iband][icol]=stat.sum(storeBuffer[bands[iband]][icol]);
1262  writeBuffer[iband][icol]=stat.sum(storeBuffer[iband][icol]);
1263  break;
1264  case(crule::minallbands):
1265  // writeBuffer[iband][icol]=stat.mymin(storeBuffer[bands[iband]][icol]);
1266  writeBuffer[iband][icol]=stat.mymin(storeBuffer[iband][icol]);
1267  break;
1268  case(crule::maxallbands):
1269  // writeBuffer[iband][icol]=stat.mymax(storeBuffer[bands[iband]][icol]);
1270  writeBuffer[iband][icol]=stat.mymax(storeBuffer[iband][icol]);
1271  break;
1272  case(crule::stdev):
1273  // writeBuffer[iband][icol]=sqrt(stat.var(storeBuffer[bands[iband]][icol]));
1274  writeBuffer[iband][icol]=sqrt(stat.var(storeBuffer[iband][icol]));
1275  break;
1276  default:
1277  break;
1278  }
1279  }
1280  catch(string error){
1281  if(verbose_opt[0])
1282  cerr << error << endl;
1283  writeBuffer[iband][icol]=dstnodata_opt[0];
1284  continue;
1285  }
1286  }
1287  try{
1288  imgWriter.writeData(writeBuffer[iband],irow,iband);
1289  }
1290  catch(string error){
1291  cerr << error << " in " << output_opt[0] << endl;
1292  throw;
1293  }
1294  }
1295  if(file_opt[0]){
1296  try{
1297  imgWriter.writeData(fileBuffer,irow,bands.size());
1298  }
1299  catch(string error){
1300  cerr << error << " in " << output_opt[0] << endl;
1301  throw;
1302  }
1303  }
1304  }
1305  progress=static_cast<float>(irow+1.0)/imgWriter.nrOfRow();
1306  pfnProgress(progress,pszMessage,pProgressArg);
1307  }
1308  if(extent_opt.size()&&(cut_opt[0]||eoption_opt.size())){
1309  extentReader.close();
1310  }
1311  for(int ifile=0;ifile<input_opt.size();++ifile){
1312  imgReader[ifile].close();
1313  }
1314  if(mask_opt.size())
1315  maskReader.close();
1316  imgWriter.close();
1317 }
1318 
void rasterizeOgr(ImgReaderOgr &ogrReader, const std::vector< double > &burnValues, const std::vector< std::string > &controlOptions=std::vector< std::string >(), const std::vector< std::string > &layernames=std::vector< std::string >())
Rasterize an OGR vector dataset using the gdal algorithm "GDALRasterizeLayers".
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
void setImageDescription(const std::string &imageDescription)
Set the image description (only for GeoTiff format: TIFFTAG_IMAGEDESCRIPTION)
Definition: ImgWriterGdal.h:55
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
CPLErr setProjectionProj4(const std::string &projection)
Set the projection for this dataset from user input (supports epsg: format) ...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96
int setNoData(const std::vector< double > nodata)
Set the no data values of this dataset using a standard template library (stl) vector as input...
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
void close(void)
Close the image.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98