| 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 int readncol=endCol-startCol+1;
684 // vector<double> readBuffer(readncol+1);
685 vector<double> readBuffer;
686 int nband=(band_opt.size())?band_opt.size() : imgReader.nrOfBand();
687 for(int iband=0;iband<nband;++iband){
688 int readBand=(band_opt.size()>iband)?band_opt[iband]:iband;
689 if(verbose_opt[0]){
690 cout << "extracting band " << readBand << endl;
691 pfnProgress(progress,pszMessage,pProgressArg);
692 }
693 double theMin=0;
694 double theMax=0;
695 if(autoscale_opt.size()){
696 try{
697 imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax);
698 }
699 catch(string errorString){
700 cout << errorString << endl;
701 }
702 if(verbose_opt[0])
703 cout << "minmax: " << theMin << ", " << theMax << endl;
704 double theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin);
705 double theOffset=autoscale_opt[0]-theScale*theMin;
706 imgReader.setScale(theScale,readBand);
707 imgReader.setOffset(theOffset,readBand);
708 }
709 else{
710 if(scale_opt.size()){
711 if(scale_opt.size()>iband)
712 imgReader.setScale(scale_opt[iband],readBand);
713 else
714 imgReader.setScale(scale_opt[0],readBand);
715 }
716 if(offset_opt.size()){
717 if(offset_opt.size()>iband)
718 imgReader.setOffset(offset_opt[iband],readBand);
719 else
720 imgReader.setOffset(offset_opt[0],readBand);
721 }
722 }
723
724 double readRow=0;
725 double readCol=0;
726 double lowerCol=0;
727 double upperCol=0;
728 for(int irow=0;irow<imgWriter.nrOfRow();++irow){
729 vector<float> lineMask;
730 double x=0;
731 double y=0;
732 //convert irow to geo
733 imgWriter.image2geo(0,irow,x,y);
734 //lookup corresponding row for irow in this file
735 imgReader.geo2image(x,y,readCol,readRow);
736 vector<double> writeBuffer;
737 if(readRow<0||readRow>=imgReader.nrOfRow()){
738 //if(readRow<0)
739 //readRow=0;
740 //else if(readRow>=imgReader.nrOfRow())
741 //readRow=imgReader.nrOfRow()-1;
742 for(int icol=0;icol<imgWriter.nrOfCol();++icol)
743 writeBuffer.push_back(nodataValue);
744 }
745 else{
746 try{
747 if(endCol<imgReader.nrOfCol()-1){
748 imgReader.readData(readBuffer,startCol,endCol+1,readRow,readBand,theResample);
749 }
750 else{
751 imgReader.readData(readBuffer,startCol,endCol,readRow,readBand,theResample);
752 }
753 // for(int icol=0;icol<ncropcol;++icol){
754 double oldRowMask=-1;//keep track of row mask to optimize number of line readings
755 for(int icol=0;icol<imgWriter.nrOfCol();++icol){
756 imgWriter.image2geo(icol,irow,x,y);
757 //lookup corresponding row for irow in this file
758 imgReader.geo2image(x,y,readCol,readRow);
759 if(readCol<0||readCol>=imgReader.nrOfCol()){
760 // if(readCol<0||readCol>=imgReader.nrOfCol()){
761 // if(readCol<0)
762 // readCol=0;
763 // else if(readCol>=imgReader.nrOfCol())
764 // readCol=imgReader.nrOfCol()-1;
765 writeBuffer.push_back(nodataValue);
766 }
767 else{
768
769 bool valid=true;
770 double geox=0;
771 double geoy=0;
772 if(mask_opt.size()){
773 //read mask
774 double colMask=0;
775 double rowMask=0;
776
777 imgWriter.image2geo(icol,irow,geox,geoy);
778 maskReader.geo2image(geox,geoy,colMask,rowMask);
779 colMask=static_cast<int>(colMask);
780 rowMask=static_cast<int>(rowMask);
781 if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
782 if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
783
784 assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
785 try{
786 maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
787 }
788 catch(string errorstring){
789 cerr << errorstring << endl;
790 exit(1);
791 }
792 catch(...){
793 cerr << "error caught" << std::endl;
794 exit(3);
795 }
796 oldRowMask=rowMask;
797 }
798 for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
799 if(lineMask[colMask]==msknodata_opt[ivalue]){
800 valid=false;
801 if(nodata_opt.size()>ivalue)
802 nodataValue=nodata_opt[ivalue];
803 }
804 }
805 }
806 }
807 if(!valid)
808 writeBuffer.push_back(nodataValue);
809 else{
810 switch(theResample){
811 case(BILINEAR):
812 lowerCol=readCol-0.5;
813 lowerCol=static_cast<int>(lowerCol);
814 upperCol=readCol+0.5;
815 upperCol=static_cast<int>(upperCol);
816 if(lowerCol<0)
817 lowerCol=0;
818 if(upperCol>=imgReader.nrOfCol())
819 upperCol=imgReader.nrOfCol()-1;
820 // writeBuffer.push_back((readCol-0.5-lowerCol)*(readBuffer[upperCol-startCol]*theScale+theOffset)+(1-readCol+0.5+lowerCol)*(readBuffer[lowerCol-startCol]*theScale+theOffset));
821 writeBuffer.push_back((readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol]);
822 break;
823 default:
824 readCol=static_cast<int>(readCol);
825 readCol-=startCol;//we only start reading from startCol
826 // writeBuffer.push_back(readBuffer[readCol]*theScale+theOffset);
827 writeBuffer.push_back(readBuffer[readCol]);
828 break;
829 }
830 }
831 }
832 }
833 }
834 catch(string errorstring){
835 cout << errorstring << endl;
836 exit(2);
837 }
838 }
839 if(writeBuffer.size()!=imgWriter.nrOfCol())
840 cout << "writeBuffer.size()=" << writeBuffer.size() << ", imgWriter.nrOfCol()=" << imgWriter.nrOfCol() << endl;
841 assert(writeBuffer.size()==imgWriter.nrOfCol());
842 try{
843 imgWriter.writeData(writeBuffer,irow,writeBand);
844 }
845 catch(string errorstring){
846 cout << errorstring << endl;
847 exit(3);
848 }
849 if(verbose_opt[0]){
850 progress=(1.0+irow);
851 progress/=imgWriter.nrOfRow();
852 pfnProgress(progress,pszMessage,pProgressArg);
853 }
854 else{
855 progress=(1.0+irow);
856 progress+=(imgWriter.nrOfRow()*writeBand);
857 progress/=imgWriter.nrOfBand()*imgWriter.nrOfRow();
858 assert(progress>=0);
859 assert(progress<=1);
860 pfnProgress(progress,pszMessage,pProgressArg);
861 }
862 }
863 ++writeBand;
864 }
865 imgReader.close();
866 }
867 if(extent_opt.size()&&(cut_opt[0]||eoption_opt.size())){
868 extentReader.close();
869 }
870 if(mask_opt.size())
871 maskReader.close();
872 imgWriter.close();
873}
Egcs
Definition: Egcs.h:27
ImgRasterGdal::getDeltaY
double getDeltaY(void) const
Get the pixel cell spacing in y.
Definition: ImgRasterGdal.h:154
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::isGeoRef
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
Definition: ImgRasterGdal.h:104
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::getProjection
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
Definition: ImgRasterGdal.cc:40
ImgRasterGdal::getUly
double getUly() const
Get the upper left corner y (georeferenced) coordinate of this dataset.
Definition: ImgRasterGdal.h:128
ImgRasterGdal::getColorTable
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
Definition: ImgRasterGdal.cc:118
ImgRasterGdal::setGeoTransform
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
Definition: ImgRasterGdal.cc:148
ImgRasterGdal::getLrx
double getLrx() const
Get the lower right corner x (georeferenced) coordinate of this dataset.
Definition: ImgRasterGdal.h:130
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::setScale
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
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::getDeltaX
double getDeltaX(void) const
Get the pixel cell spacing in x.
Definition: ImgRasterGdal.h:152
ImgRasterGdal::setOffset
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
ImgRasterGdal::getInterleave
std::string getInterleave() const
Get the band coding (interleave)
Definition: ImgRasterGdal.cc:293
ImgRasterGdal::getLry
double getLry() const
Get the lower right corner y (georeferenced) coordinate of this dataset.
Definition: ImgRasterGdal.h:132
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
ImgRasterGdal::getUlx
double getUlx() const
Get the upper left corner x (georeferenced) coordinate of this dataset.
Definition: ImgRasterGdal.h:126
ImgRasterGdal::getBoundingBox
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
Definition: ImgRasterGdal.cc:333
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::getMinMax
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...
Definition: ImgReaderGdal.cc:161
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