pktools  2.6.7
Processing Kernel for geospatial data
pkcrop.cc
1 /**********************************************************************
2 pkcrop.cc: perform raster data operations on image such as crop, extract and stack bands
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 <assert.h>
21 #include <cstdlib>
22 #include <string>
23 #include <list>
24 #include <iostream>
25 #include <algorithm>
26 #include "imageclasses/ImgWriterGdal.h"
27 #include "imageclasses/ImgReaderGdal.h"
28 #include "imageclasses/ImgReaderOgr.h"
29 #include "base/Optionpk.h"
30 #include "algorithms/Egcs.h"
31 
32 /******************************************************************************/
99 using namespace std;
100 
101 int main(int argc, char *argv[])
102 {
103  Optionpk<string> input_opt("i", "input", "Input image file(s). If input contains multiple images, a multi-band output is created");
104  Optionpk<string> output_opt("o", "output", "Output image file");
105  Optionpk<string> projection_opt("a_srs", "a_srs", "Override the projection for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid");
106  //todo: support layer names
107  Optionpk<string> extent_opt("e", "extent", "get boundary from extent from polygons in vector file");
108  Optionpk<bool> cut_opt("cut", "crop_to_cutline", "Crop the extent of the target dataset to the extent of the cutline.",false);
109  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");
110  Optionpk<string> mask_opt("m", "mask", "Use the the specified file as a validity mask (0 is nodata).");
111  Optionpk<float> msknodata_opt("msknodata", "msknodata", "Mask value not to consider for crop.", 0);
112  Optionpk<short> mskband_opt("mskband", "mskband", "Mask band to read (0 indexed)", 0);
113  Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box", 0.0);
114  Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box", 0.0);
115  Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box", 0.0);
116  Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box", 0.0);
117  Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (empty: keep original resolution)");
118  Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (empty: keep original resolution)");
119  Optionpk<double> cx_opt("x", "x", "x-coordinate of image center to crop (in meter)");
120  Optionpk<double> cy_opt("y", "y", "y-coordinate of image center to crop (in meter)");
121  Optionpk<double> nx_opt("nx", "nx", "image size in x to crop (in meter)");
122  Optionpk<double> ny_opt("ny", "ny", "image size in y to crop (in meter)");
123  Optionpk<int> ns_opt("ns", "ns", "number of samples to crop (in pixels)");
124  Optionpk<int> nl_opt("nl", "nl", "number of lines to crop (in pixels)");
125  Optionpk<unsigned short> band_opt("b", "band", "band index to crop (leave empty to retain all bands)");
126  Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
127  Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
128  Optionpk<double> autoscale_opt("as", "autoscale", "scale output to min and max, e.g., --autoscale 0 --autoscale 255");
129  Optionpk<double> scale_opt("scale", "scale", "output=scale*input+offset");
130  Optionpk<double> offset_opt("offset", "offset", "output=scale*input+offset");
131  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","");
132  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
133  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
134  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
135  Optionpk<double> nodata_opt("nodata", "nodata", "Nodata value to put in image if out of bounds.");
136  Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation).", "near");
137  Optionpk<string> description_opt("d", "description", "Set image description");
138  Optionpk<bool> align_opt("align", "align", "Align output bounding box to input image",false);
139  Optionpk<short> verbose_opt("v", "verbose", "verbose", 0,2);
140 
141  extent_opt.setHide(1);
142  cut_opt.setHide(1);
143  eoption_opt.setHide(1);
144  bstart_opt.setHide(1);
145  bend_opt.setHide(1);
146  mask_opt.setHide(1);
147  msknodata_opt.setHide(1);
148  mskband_opt.setHide(1);
149  option_opt.setHide(1);
150  cx_opt.setHide(1);
151  cy_opt.setHide(1);
152  nx_opt.setHide(1);
153  ny_opt.setHide(1);
154  ns_opt.setHide(1);
155  nl_opt.setHide(1);
156  scale_opt.setHide(1);
157  offset_opt.setHide(1);
158  nodata_opt.setHide(1);
159  description_opt.setHide(1);
160 
161  bool doProcess;//stop process when program was invoked with help option (-h --help)
162  try{
163  doProcess=input_opt.retrieveOption(argc,argv);
164  output_opt.retrieveOption(argc,argv);
165  projection_opt.retrieveOption(argc,argv);
166  ulx_opt.retrieveOption(argc,argv);
167  uly_opt.retrieveOption(argc,argv);
168  lrx_opt.retrieveOption(argc,argv);
169  lry_opt.retrieveOption(argc,argv);
170  band_opt.retrieveOption(argc,argv);
171  bstart_opt.retrieveOption(argc,argv);
172  bend_opt.retrieveOption(argc,argv);
173  autoscale_opt.retrieveOption(argc,argv);
174  otype_opt.retrieveOption(argc,argv);
175  oformat_opt.retrieveOption(argc,argv);
176  colorTable_opt.retrieveOption(argc,argv);
177  dx_opt.retrieveOption(argc,argv);
178  dy_opt.retrieveOption(argc,argv);
179  resample_opt.retrieveOption(argc,argv);
180  extent_opt.retrieveOption(argc,argv);
181  cut_opt.retrieveOption(argc,argv);
182  eoption_opt.retrieveOption(argc,argv);
183  mask_opt.retrieveOption(argc,argv);
184  msknodata_opt.retrieveOption(argc,argv);
185  mskband_opt.retrieveOption(argc,argv);
186  option_opt.retrieveOption(argc,argv);
187  cx_opt.retrieveOption(argc,argv);
188  cy_opt.retrieveOption(argc,argv);
189  nx_opt.retrieveOption(argc,argv);
190  ny_opt.retrieveOption(argc,argv);
191  ns_opt.retrieveOption(argc,argv);
192  nl_opt.retrieveOption(argc,argv);
193  scale_opt.retrieveOption(argc,argv);
194  offset_opt.retrieveOption(argc,argv);
195  nodata_opt.retrieveOption(argc,argv);
196  description_opt.retrieveOption(argc,argv);
197  align_opt.retrieveOption(argc,argv);
198  verbose_opt.retrieveOption(argc,argv);
199  }
200  catch(string predefinedString){
201  std::cout << predefinedString << std::endl;
202  exit(0);
203  }
204  if(verbose_opt[0])
205  cout << setprecision(12) << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
206 
207  if(!doProcess){
208  cout << endl;
209  cout << "Usage: pkcrop -i input -o output" << endl;
210  cout << endl;
211  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
212  exit(0);//help was invoked, stop processing
213  }
214  if(input_opt.empty()){
215  std::cerr << "No input file provided (use option -i). Use --help for help information" << std::endl;
216  exit(0);
217  }
218  if(output_opt.empty()){
219  std::cerr << "No output file provided (use option -o). Use --help for help information" << std::endl;
220  exit(0);
221  }
222 
223  float nodataValue=nodata_opt.size()? nodata_opt[0] : 0;
224  RESAMPLE theResample;
225  if(resample_opt[0]=="near"){
226  theResample=NEAR;
227  if(verbose_opt[0])
228  cout << "resampling: nearest neighbor" << endl;
229  }
230  else if(resample_opt[0]=="bilinear"){
231  theResample=BILINEAR;
232  if(verbose_opt[0])
233  cout << "resampling: bilinear interpolation" << endl;
234  }
235  else{
236  std::cout << "Error: resampling method " << resample_opt[0] << " not supported" << std::endl;
237  exit(1);
238  }
239 
240  const char* pszMessage;
241  void* pProgressArg=NULL;
242  GDALProgressFunc pfnProgress=GDALTermProgress;
243  double progress=0;
244  pfnProgress(progress,pszMessage,pProgressArg);
245  ImgReaderGdal imgReader;
246  ImgWriterGdal imgWriter;
247  //open input images to extract number of bands and spatial resolution
248  int ncropband=0;//total number of bands to write
249  double dx=0;
250  double dy=0;
251  if(dx_opt.size())
252  dx=dx_opt[0];
253  if(dy_opt.size())
254  dy=dy_opt[0];
255 
256  //convert start and end band options to vector of band indexes
257  try{
258  if(bstart_opt.size()){
259  if(bend_opt.size()!=bstart_opt.size()){
260  string errorstring="Error: options for start and end band indexes must be provided as pairs, missing end band";
261  throw(errorstring);
262  }
263  band_opt.clear();
264  for(int ipair=0;ipair<bstart_opt.size();++ipair){
265  if(bend_opt[ipair]<=bstart_opt[ipair]){
266  string errorstring="Error: index for end band must be smaller then start band";
267  throw(errorstring);
268  }
269  for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
270  band_opt.push_back(iband);
271  }
272  }
273  }
274  catch(string error){
275  cerr << error << std::endl;
276  exit(1);
277  }
278 
279  bool isGeoRef=false;
280  string projectionString;
281  for(int iimg=0;iimg<input_opt.size();++iimg){
282  try{
283  imgReader.open(input_opt[iimg],GA_ReadOnly);
284  }
285  catch(string error){
286  cerr << "Error: could not open file " << input_opt[iimg] << ": " << error << std::endl;
287  exit(1);
288  }
289  if(!isGeoRef)
290  isGeoRef=imgReader.isGeoRef();
291  if(imgReader.isGeoRef()&&projection_opt.empty())
292  projectionString=imgReader.getProjection();
293  if(dx_opt.empty()){
294  if(!iimg||imgReader.getDeltaX()<dx)
295  dx=imgReader.getDeltaX();
296  }
297 
298  if(dy_opt.empty()){
299  if(!iimg||imgReader.getDeltaY()<dy)
300  dy=imgReader.getDeltaY();
301  }
302  if(band_opt.size())
303  ncropband+=band_opt.size();
304  else
305  ncropband+=imgReader.nrOfBand();
306  imgReader.close();
307  }
308 
309  GDALDataType theType=GDT_Unknown;
310  if(verbose_opt[0])
311  cout << "possible output data types: ";
312  for(int iType = 0; iType < GDT_TypeCount; ++iType){
313  if(verbose_opt[0])
314  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
315  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
316  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
317  otype_opt[0].c_str()))
318  theType=(GDALDataType) iType;
319  }
320  if(verbose_opt[0]){
321  cout << endl;
322  if(theType==GDT_Unknown)
323  cout << "Unknown output pixel type: " << otype_opt[0] << endl;
324  else
325  cout << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
326  }
327  //bounding box of cropped image
328  double cropulx=ulx_opt[0];
329  double cropuly=uly_opt[0];
330  double croplrx=lrx_opt[0];
331  double croplry=lry_opt[0];
332  //get bounding box from extentReader if defined
333  ImgReaderOgr extentReader;
334 
335  if(extent_opt.size()){
336  double e_ulx;
337  double e_uly;
338  double e_lrx;
339  double e_lry;
340  for(int iextent=0;iextent<extent_opt.size();++iextent){
341  try{
342  extentReader.open(extent_opt[iextent]);
343  if(!(extentReader.getExtent(e_ulx,e_uly,e_lrx,e_lry))){
344  ostringstream os;
345  os << "Error: could not get extent from " << extent_opt[0] << endl;
346  throw(os.str());
347  }
348  }
349  catch(string error){
350  cerr << error << std::endl;
351  exit(1);
352  }
353  if(!iextent){
354  ulx_opt[0]=e_ulx;
355  uly_opt[0]=e_uly;
356  lrx_opt[0]=e_lrx;
357  lry_opt[0]=e_lry;
358  }
359  else{
360  if(e_ulx<ulx_opt[0])
361  ulx_opt[0]=e_ulx;
362  if(e_uly>uly_opt[0])
363  uly_opt[0]=e_uly;
364  if(e_lrx>lrx_opt[0])
365  lrx_opt[0]=e_lrx;
366  if(e_lry<lry_opt[0])
367  lry_opt[0]=e_lry;
368  }
369  extentReader.close();
370  }
371  if(croplrx>cropulx&&cropulx>ulx_opt[0])
372  ulx_opt[0]=cropulx;
373  if(croplrx>cropulx&&croplrx<lrx_opt[0])
374  lrx_opt[0]=croplrx;
375  if(cropuly>croplry&&cropuly<uly_opt[0])
376  uly_opt[0]=cropuly;
377  if(croplry<cropuly&&croplry>lry_opt[0])
378  lry_opt[0]=croplry;
379  if(cut_opt.size()||eoption_opt.size())
380  extentReader.open(extent_opt[0]);
381  }
382  else if(cx_opt.size()&&cy_opt.size()&&nx_opt.size()&&ny_opt.size()){
383  ulx_opt[0]=cx_opt[0]-nx_opt[0]/2.0;
384  uly_opt[0]=(isGeoRef) ? cy_opt[0]+ny_opt[0]/2.0 : cy_opt[0]-ny_opt[0]/2.0;
385  lrx_opt[0]=cx_opt[0]+nx_opt[0]/2.0;
386  lry_opt[0]=(isGeoRef) ? cy_opt[0]-ny_opt[0]/2.0 : cy_opt[0]+ny_opt[0]/2.0;
387  // if(cropulx<ulx_opt[0])
388  // cropulx=ulx_opt[0];
389  // if(cropuly>uly_opt[0])
390  // cropuly=uly_opt[0];
391  // if(croplrx>lrx_opt[0])
392  // croplrx=lrx_opt[0];
393  // if(croplry<lry_opt[0])
394  // croplry=lry_opt[0];
395  }
396  else if(cx_opt.size()&&cy_opt.size()&&ns_opt.size()&&nl_opt.size()){
397  ulx_opt[0]=cx_opt[0]-ns_opt[0]*dx/2.0;
398  uly_opt[0]=(isGeoRef) ? cy_opt[0]+nl_opt[0]*dy/2.0 : cy_opt[0]-nl_opt[0]*dy/2.0;
399  lrx_opt[0]=cx_opt[0]+ns_opt[0]*dx/2.0;
400  lry_opt[0]=(isGeoRef) ? cy_opt[0]-nl_opt[0]*dy/2.0 : cy_opt[0]+nl_opt[0]*dy/2.0;
401  // if(cropulx<ulx_opt[0])
402  // cropulx=ulx_opt[0];
403  // if(cropuly>uly_opt[0])
404  // cropuly=uly_opt[0];
405  // if(croplrx>lrx_opt[0])
406  // croplrx=lrx_opt[0];
407  // if(croplry<lry_opt[0])
408  // croplry=lry_opt[0];
409  }
410 
411  if(verbose_opt[0])
412  cout << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
413 
414  int ncropcol=0;
415  int ncroprow=0;
416 
417  ImgWriterGdal maskWriter;
418  if(extent_opt.size()&&(cut_opt[0]||eoption_opt.size())){
419  if(mask_opt.size()){
420  string errorString="Error: can only either mask or extent extent with cutline, not both";
421  throw(errorString);
422  }
423  try{
424  ncropcol=abs(static_cast<int>(ceil((lrx_opt[0]-ulx_opt[0])/dx)));
425  ncroprow=abs(static_cast<int>(ceil((uly_opt[0]-lry_opt[0])/dy)));
426  //todo: produce unique name
427  maskWriter.open("/vsimem/mask.tif",ncropcol,ncroprow,1,GDT_Float32,"GTiff",option_opt);
428  double gt[6];
429  gt[0]=ulx_opt[0];
430  gt[1]=dx;
431  gt[2]=0;
432  gt[3]=uly_opt[0];
433  gt[4]=0;
434  gt[5]=-dy;
435  maskWriter.setGeoTransform(gt);
436  if(projection_opt.size())
437  maskWriter.setProjectionProj4(projection_opt[0]);
438  else if(projectionString.size())
439  maskWriter.setProjection(projectionString);
440 
441  vector<double> burnValues(1,1);//burn value is 1 (single band)
442  maskWriter.rasterizeOgr(extentReader,burnValues,eoption_opt);
443  maskWriter.close();
444  }
445  catch(string error){
446  cerr << error << std::endl;
447  exit(2);
448  }
449  mask_opt.clear();
450  mask_opt.push_back("/vsimem/mask.tif");
451  }
452  ImgReaderGdal maskReader;
453  if(mask_opt.size()==1){
454  try{
455  //there is only a single mask
456  maskReader.open(mask_opt[0],GA_ReadOnly);
457  if(mskband_opt[0]>=maskReader.nrOfBand()){
458  string errorString="Error: illegal mask band";
459  throw(errorString);
460  }
461  }
462  catch(string error){
463  cerr << error << std::endl;
464  exit(2);
465  }
466  }
467 
468  //determine number of output bands
469  int writeBand=0;//write band
470 
471  if(scale_opt.size()){
472  while(scale_opt.size()<band_opt.size())
473  scale_opt.push_back(scale_opt[0]);
474  }
475  if(offset_opt.size()){
476  while(offset_opt.size()<band_opt.size())
477  offset_opt.push_back(offset_opt[0]);
478  }
479  if(autoscale_opt.size()){
480  assert(autoscale_opt.size()%2==0);
481  // while(autoscale_opt.size()<band_opt.size()*2){
482  // autoscale_opt.push_back(autoscale_opt[0]);
483  // autoscale_opt.push_back(autoscale_opt[1]);
484  // }
485  }
486 
487  for(int iimg=0;iimg<input_opt.size();++iimg){
488  if(verbose_opt[0])
489  cout << "opening image " << input_opt[iimg] << endl;
490  try{
491  imgReader.open(input_opt[iimg],GA_ReadOnly);
492  }
493  catch(string error){
494  cerr << error << std::endl;
495  exit(2);
496  }
497  //if output type not set, get type from input image
498  if(theType==GDT_Unknown){
499  theType=imgReader.getDataType();
500  if(verbose_opt[0])
501  cout << "Using data type from input image: " << GDALGetDataTypeName(theType) << endl;
502  }
503  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
504  string theInterleave="INTERLEAVE=";
505  theInterleave+=imgReader.getInterleave();
506  option_opt.push_back(theInterleave);
507  }
508  int nrow=imgReader.nrOfRow();
509  int ncol=imgReader.nrOfCol();
510  // if(!dx||!dy){
511  // dx=imgReader.getDeltaX();
512  // dy=imgReader.getDeltaY();
513  // }
514  if(verbose_opt[0])
515  cout << "size of " << input_opt[iimg] << ": " << ncol << " cols, "<< nrow << " rows" << endl;
516  double uli,ulj,lri,lrj;//image coordinates
517  bool forceEUgrid=false;
518  if(projection_opt.size())
519  forceEUgrid=(!(projection_opt[0].compare("EPSG:3035"))||!(projection_opt[0].compare("EPSG:3035"))||projection_opt[0].find("ETRS-LAEA")!=string::npos);
520  if(ulx_opt[0]>=lrx_opt[0]){//default bounding box: no cropping
521  uli=0;
522  lri=imgReader.nrOfCol()-1;
523  ulj=0;
524  lrj=imgReader.nrOfRow()-1;
525  ncropcol=imgReader.nrOfCol();
526  ncroprow=imgReader.nrOfRow();
527  imgReader.getBoundingBox(cropulx,cropuly,croplrx,croplry);
528  double magicX=1,magicY=1;
529  // imgReader.getMagicPixel(magicX,magicY);
530  if(forceEUgrid){
531  //force to LAEA grid
532  Egcs egcs;
533  egcs.setLevel(egcs.res2level(dx));
534  egcs.force2grid(cropulx,cropuly,croplrx,croplry);
535  imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
536  imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
537  }
538  imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
539  imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
540  //test
541  ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx)));
542  ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy)));
543  }
544  else{
545  double magicX=1,magicY=1;
546  // imgReader.getMagicPixel(magicX,magicY);
547  cropulx=ulx_opt[0];
548  cropuly=uly_opt[0];
549  croplrx=lrx_opt[0];
550  croplry=lry_opt[0];
551  if(forceEUgrid){
552  //force to LAEA grid
553  Egcs egcs;
554  egcs.setLevel(egcs.res2level(dx));
555  egcs.force2grid(cropulx,cropuly,croplrx,croplry);
556  }
557  else if(align_opt[0]){
558  if(cropulx>imgReader.getUlx())
559  cropulx-=fmod(cropulx-imgReader.getUlx(),dx);
560  else if(cropulx<imgReader.getUlx())
561  cropulx+=fmod(imgReader.getUlx()-cropulx,dx)-dx;
562  if(croplrx<imgReader.getLrx())
563  croplrx+=fmod(imgReader.getLrx()-croplrx,dx);
564  else if(croplrx>imgReader.getLrx())
565  croplrx-=fmod(croplrx-imgReader.getLrx(),dx)+dx;
566  if(croplry>imgReader.getLry())
567  croplry-=fmod(croplry-imgReader.getLry(),dy);
568  else if(croplry<imgReader.getLry())
569  croplry+=fmod(imgReader.getLry()-croplry,dy)-dy;
570  if(cropuly<imgReader.getUly())
571  cropuly+=fmod(imgReader.getUly()-cropuly,dy);
572  else if(cropuly>imgReader.getUly())
573  cropuly-=fmod(cropuly-imgReader.getUly(),dy)+dy;
574  }
575  imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
576  imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
577 
578  ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx)));
579  ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy)));
580  uli=floor(uli);
581  ulj=floor(ulj);
582  lri=floor(lri);
583  lrj=floor(lrj);
584  }
585 
586  double dcropcol=0;
587  double dcroprow=0;
588  double deltaX=imgReader.getDeltaX();
589  double deltaY=imgReader.getDeltaY();
590  dcropcol=(lri-uli+1)/(dx/deltaX);
591  dcroprow=(lrj-ulj+1)/(dy/deltaY);
592  if(!imgWriter.nrOfBand()){//not opened yet
593  if(verbose_opt[0]){
594  cout << "cropulx: " << cropulx << endl;
595  cout << "cropuly: " << cropuly << endl;
596  cout << "croplrx: " << croplrx << endl;
597  cout << "croplry: " << croplry << endl;
598  cout << "ncropcol: " << ncropcol << endl;
599  cout << "ncroprow: " << ncroprow << endl;
600  cout << "cropulx+ncropcol*dx: " << cropulx+ncropcol*dx << endl;
601  cout << "cropuly-ncroprow*dy: " << cropuly-ncroprow*dy << endl;
602  cout << "upper left column of input image: " << uli << endl;
603  cout << "upper left row of input image: " << ulj << endl;
604  cout << "lower right column of input image: " << lri << endl;
605  cout << "lower right row of input image: " << lrj << endl;
606  cout << "new number of cols: " << ncropcol << endl;
607  cout << "new number of rows: " << ncroprow << endl;
608  cout << "new number of bands: " << ncropband << endl;
609  }
610  // string theCompression;
611  // if(compress_opt[0]!="")//default
612  // theCompression=compress_opt[0];
613  // else
614  // theCompression=imgReader.getCompression();
615  // string theInterleave;
616  // if(interleave_opt[0]!="")//default
617  // theInterleave=interleave_opt[0];
618  // else
619  // theInterleave=imgReader.getInterleave();
620  string imageType;//=imgReader.getImageType();
621  if(oformat_opt.size())//default
622  imageType=oformat_opt[0];
623  try{
624  imgWriter.open(output_opt[0],ncropcol,ncroprow,ncropband,theType,imageType,option_opt);
625  if(nodata_opt.size()){
626  imgWriter.setNoData(nodata_opt);
627  // for(int iband=0;iband<ncropband;++iband)
628  // imgWriter.GDALSetNoDataValue(nodata_opt[0],iband);
629  }
630  }
631  catch(string errorstring){
632  cout << errorstring << endl;
633  exit(4);
634  }
635  if(description_opt.size())
636  imgWriter.setImageDescription(description_opt[0]);
637  double gt[6];
638  gt[0]=cropulx;
639  gt[1]=dx;
640  gt[2]=0;
641  gt[3]=cropuly;
642  gt[4]=0;
643  gt[5]=(imgReader.isGeoRef())? -dy : dy;
644  imgWriter.setGeoTransform(gt);
645  if(projection_opt.size()){
646  if(verbose_opt[0])
647  cout << "projection: " << projection_opt[0] << endl;
648  imgWriter.setProjectionProj4(projection_opt[0]);
649  }
650  else
651  imgWriter.setProjection(imgReader.getProjection());
652  if(imgWriter.getDataType()==GDT_Byte){
653  if(colorTable_opt.size()){
654  if(colorTable_opt[0]!="none")
655  imgWriter.setColorTable(colorTable_opt[0]);
656  }
657  else if (imgReader.getColorTable()!=NULL)//copy colorTable from input image
658  imgWriter.setColorTable(imgReader.getColorTable());
659  }
660  }
661 
662  double startCol=uli;
663  double endCol=lri;
664  if(uli<0)
665  startCol=0;
666  else if(uli>=imgReader.nrOfCol())
667  startCol=imgReader.nrOfCol()-1;
668  if(lri<0)
669  endCol=0;
670  else if(lri>=imgReader.nrOfCol())
671  endCol=imgReader.nrOfCol()-1;
672  double startRow=ulj;
673  double endRow=lrj;
674  if(ulj<0)
675  startRow=0;
676  else if(ulj>=imgReader.nrOfRow())
677  startRow=imgReader.nrOfRow()-1;
678  if(lrj<0)
679  endRow=0;
680  else if(lrj>=imgReader.nrOfRow())
681  endRow=imgReader.nrOfRow()-1;
682 
683 
684 
685  int readncol=endCol-startCol+1;
686  // vector<double> readBuffer(readncol+1);
687  vector<double> readBuffer;
688  int nband=(band_opt.size())?band_opt.size() : imgReader.nrOfBand();
689  for(int iband=0;iband<nband;++iband){
690  int readBand=(band_opt.size()>iband)?band_opt[iband]:iband;
691  if(verbose_opt[0]){
692  cout << "extracting band " << readBand << endl;
693  pfnProgress(progress,pszMessage,pProgressArg);
694  }
695  double theMin=0;
696  double theMax=0;
697  if(autoscale_opt.size()){
698  try{
699  imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax);
700  }
701  catch(string errorString){
702  cout << errorString << endl;
703  }
704  if(verbose_opt[0])
705  cout << "minmax: " << theMin << ", " << theMax << endl;
706  double theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin);
707  double theOffset=autoscale_opt[0]-theScale*theMin;
708  imgReader.setScale(theScale,readBand);
709  imgReader.setOffset(theOffset,readBand);
710  }
711  else{
712  if(scale_opt.size()){
713  if(scale_opt.size()>iband)
714  imgReader.setScale(scale_opt[iband],readBand);
715  else
716  imgReader.setScale(scale_opt[0],readBand);
717  }
718  if(offset_opt.size()){
719  if(offset_opt.size()>iband)
720  imgReader.setOffset(offset_opt[iband],readBand);
721  else
722  imgReader.setOffset(offset_opt[0],readBand);
723  }
724  }
725 
726  double readRow=0;
727  double readCol=0;
728  double lowerCol=0;
729  double upperCol=0;
730  for(int irow=0;irow<imgWriter.nrOfRow();++irow){
731  vector<float> lineMask;
732  double x=0;
733  double y=0;
734  //convert irow to geo
735  imgWriter.image2geo(0,irow,x,y);
736  //lookup corresponding row for irow in this file
737  imgReader.geo2image(x,y,readCol,readRow);
738  vector<double> writeBuffer;
739  if(readRow<0||readRow>=imgReader.nrOfRow()){
740  //if(readRow<0)
741  //readRow=0;
742  //else if(readRow>=imgReader.nrOfRow())
743  //readRow=imgReader.nrOfRow()-1;
744  for(int icol=0;icol<imgWriter.nrOfCol();++icol)
745  writeBuffer.push_back(nodataValue);
746  }
747  else{
748  try{
749  if(endCol<imgReader.nrOfCol()-1){
750  imgReader.readData(readBuffer,startCol,endCol+1,readRow,readBand,theResample);
751  }
752  else{
753  imgReader.readData(readBuffer,startCol,endCol,readRow,readBand,theResample);
754  }
755  // for(int icol=0;icol<ncropcol;++icol){
756  double oldRowMask=-1;//keep track of row mask to optimize number of line readings
757  for(int icol=0;icol<imgWriter.nrOfCol();++icol){
758  imgWriter.image2geo(icol,irow,x,y);
759  //lookup corresponding row for irow in this file
760  imgReader.geo2image(x,y,readCol,readRow);
761  if(readCol<0||readCol>=imgReader.nrOfCol()){
762  // if(readCol<0||readCol>=imgReader.nrOfCol()){
763  // if(readCol<0)
764  // readCol=0;
765  // else if(readCol>=imgReader.nrOfCol())
766  // readCol=imgReader.nrOfCol()-1;
767  writeBuffer.push_back(nodataValue);
768  }
769  else{
770 
771  bool valid=true;
772  double geox=0;
773  double geoy=0;
774  if(mask_opt.size()){
775  //read mask
776  double colMask=0;
777  double rowMask=0;
778 
779  imgWriter.image2geo(icol,irow,geox,geoy);
780  maskReader.geo2image(geox,geoy,colMask,rowMask);
781  colMask=static_cast<int>(colMask);
782  rowMask=static_cast<int>(rowMask);
783  if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
784  if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
785 
786  assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
787  try{
788  maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
789  }
790  catch(string errorstring){
791  cerr << errorstring << endl;
792  exit(1);
793  }
794  catch(...){
795  cerr << "error caught" << std::endl;
796  exit(3);
797  }
798  oldRowMask=rowMask;
799  }
800  for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
801  if(lineMask[colMask]==msknodata_opt[ivalue]){
802  valid=false;
803  if(nodata_opt.size()>ivalue)
804  nodataValue=nodata_opt[ivalue];
805  }
806  }
807  }
808  }
809  if(!valid)
810  writeBuffer.push_back(nodataValue);
811  else{
812  switch(theResample){
813  case(BILINEAR):
814  lowerCol=readCol-0.5;
815  lowerCol=static_cast<int>(lowerCol);
816  upperCol=readCol+0.5;
817  upperCol=static_cast<int>(upperCol);
818  if(lowerCol<0)
819  lowerCol=0;
820  if(upperCol>=imgReader.nrOfCol())
821  upperCol=imgReader.nrOfCol()-1;
822  // writeBuffer.push_back((readCol-0.5-lowerCol)*(readBuffer[upperCol-startCol]*theScale+theOffset)+(1-readCol+0.5+lowerCol)*(readBuffer[lowerCol-startCol]*theScale+theOffset));
823  writeBuffer.push_back((readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol]);
824  break;
825  default:
826  readCol=static_cast<int>(readCol);
827  readCol-=startCol;//we only start reading from startCol
828  // writeBuffer.push_back(readBuffer[readCol]*theScale+theOffset);
829  writeBuffer.push_back(readBuffer[readCol]);
830  break;
831  }
832  }
833  }
834  }
835  }
836  catch(string errorstring){
837  cout << errorstring << endl;
838  exit(2);
839  }
840  }
841  if(writeBuffer.size()!=imgWriter.nrOfCol())
842  cout << "writeBuffer.size()=" << writeBuffer.size() << ", imgWriter.nrOfCol()=" << imgWriter.nrOfCol() << endl;
843  assert(writeBuffer.size()==imgWriter.nrOfCol());
844  try{
845  imgWriter.writeData(writeBuffer,irow,writeBand);
846  }
847  catch(string errorstring){
848  cout << errorstring << endl;
849  exit(3);
850  }
851  if(verbose_opt[0]){
852  progress=(1.0+irow);
853  progress/=imgWriter.nrOfRow();
854  pfnProgress(progress,pszMessage,pProgressArg);
855  }
856  else{
857  progress=(1.0+irow);
858  progress+=(imgWriter.nrOfRow()*writeBand);
859  progress/=imgWriter.nrOfBand()*imgWriter.nrOfRow();
860  assert(progress>=0);
861  assert(progress<=1);
862  pfnProgress(progress,pszMessage,pProgressArg);
863  }
864  }
865  ++writeBand;
866  }
867  imgReader.close();
868  }
869  if(extent_opt.size()&&(cut_opt[0]||eoption_opt.size())){
870  extentReader.close();
871  }
872  if(mask_opt.size())
873  maskReader.close();
874  imgWriter.close();
875 }
double getLry() const
Get the lower right corner y (georeferenced) coordinate of this dataset.
void setOffset(double theOffset, int band=0)
Set offset for a specific band when writing the raster data values. The scaling and offset are applie...
Definition: ImgRasterGdal.h:84
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) ...
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
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.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
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)
double getUly() const
Get the upper left corner y (georeferenced) coordinate of this dataset.
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
double getDeltaY(void) const
Get the pixel cell spacing in y.
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
double getLrx() const
Get the lower right corner x (georeferenced) coordinate of 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 getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double &minValue, double &maxValue)
Get the minimum and maximum cell values for a specific band in a region of interest defined by startC...
double getUlx() const
Get the upper left corner x (georeferenced) coordinate of this dataset.
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...
double getDeltaX(void) const
Get the pixel cell spacing in x.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
Definition: Egcs.h:26
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
void setScale(double theScale, int band=0)
Set scale for a specific band when writing the raster data values. The scaling and offset are applied...
Definition: ImgRasterGdal.h:75
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98