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