| pktools 2.6.7 Processing Kernel for geospatial data

pkdumpimg.cc
1/**********************************************************************
2 pkdumpimg.cc: program to dump image content to ascii or std out
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 <string>
21#include <fstream>
22#include <vector>
23#include <iostream>
24#include <assert.h>
25#include "base/Optionpk.h"
26#include "imageclasses/ImgReaderOgr.h"
27#include "imageclasses/ImgWriterGdal.h"
28// #include "imageclasses/ImgWriterOgr.h"
29
30#ifdef HAVE_CONFIG_H
31#include <config.h>
32#endif
33
34/******************************************************************************/
84 using namespace std;
85
86 int main(int argc, char *argv[])
87{
88 Optionpk<std::string> input_opt("i","input","input image file","");
89 Optionpk<string> output_opt("o", "output", "Output ascii file (Default is empty: use stdout", "");
90 Optionpk<string> oformat_opt("of", "oformat", "Output format (matrix form or list (x,y,z) form. Default is matrix form", "matrix");
91 Optionpk<int> band_opt("b", "band", "band index to crop");
92 Optionpk<string> extent_opt("e", "extent", "get boundary from extent from polygons in vector file", "");
93 Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box (in geocoordinates if georef is true)",0.0);
94 Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box (in geocoordinates if georef is true)",0.0);
95 Optionpk<double> lrx_opt("lrx", "lrx", "Lower left x value bounding box (in geocoordinates if georef is true)",0.0);
96 Optionpk<double> lry_opt("lry", "lry", "Lower left y value bounding box (in geocoordinates if georef is true)",0.0);
97 Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (0.0: keep original resolution)",0.0);
98 Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (0.0: keep original resolution)",0.0);
99 Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
100 Optionpk<short> dstnodata_opt("dstnodata", "dstnodata", "nodata value for output if out of bounds.", 0);
101 Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "set no data value(s) for input image");
102 Optionpk<short> verbose_opt("v", "verbose", "verbose (Default: 0)", 0,2);
103
104 dx_opt.setHide(1);
105 dy_opt.setHide(1);
106 resample_opt.setHide(1);
107 srcnodata_opt.setHide(1);
108 dstnodata_opt.setHide(1);
109
110 bool doProcess;//stop process when program was invoked with help option (-h --help)
111 try{
112 doProcess=input_opt.retrieveOption(argc,argv);
113 output_opt.retrieveOption(argc,argv);
114 oformat_opt.retrieveOption(argc,argv);
115 band_opt.retrieveOption(argc,argv);
116 extent_opt.retrieveOption(argc,argv);
117 ulx_opt.retrieveOption(argc,argv);
118 uly_opt.retrieveOption(argc,argv);
119 lrx_opt.retrieveOption(argc,argv);
120 lry_opt.retrieveOption(argc,argv);
121 dx_opt.retrieveOption(argc,argv);
122 dy_opt.retrieveOption(argc,argv);
123 resample_opt.retrieveOption(argc,argv);
124 srcnodata_opt.retrieveOption(argc,argv);
125 dstnodata_opt.retrieveOption(argc,argv);
126 verbose_opt.retrieveOption(argc,argv);
127 }
128 catch(string predefinedString){
129 std::cout << predefinedString << std::endl;
130 exit(0);
131 }
132 if(!doProcess){
133 cout << endl;
134 cout << "Usage: pkdumpimg -i input.txt [-o output]" << endl;
135 cout << endl;
136 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
137 exit(0);//help was invoked, stop processing
138 }
139
140 ofstream outputStream;
141 if(!output_opt.empty())
142 outputStream.open(output_opt[0].c_str());
143
144 RESAMPLE theResample;
145 if(resample_opt[0]=="near"){
146 theResample=NEAR;
147 if(verbose_opt[0])
148 cout << "resampling: nearest neighbour" << endl;
149 }
150 else if(resample_opt[0]=="bilinear"){
151 theResample=BILINEAR;
152 if(verbose_opt[0])
153 cout << "resampling: bilinear interpolation" << endl;
154 }
155 else{
156 std::cout << "Error: resampling method " << resample_opt[0] << " not supported" << std::endl;
157 exit(1);
158 }
159
160 // ImgWriterGdal imgWriter;
161 GDALDataType theType;
162
163 if(input_opt.empty()){
164 std::cerr << "No input file provided (use option -i). Use --help for help information" << std::endl;
165 exit(0);
166 }
167
168 ImgReaderGdal imgReader(input_opt[0]);
169 for(int inodata=0;inodata<srcnodata_opt.size();++inodata)
170 imgReader.pushNoDataValue(srcnodata_opt[inodata]);
171
172 // ImgWriterGdal virtualWriter;//only for coordinate conversion (no output file defined)
173
174 int nband=imgReader.nrOfBand();
175 //get number of lines
176 int nrow=imgReader.nrOfRow();
177 int ncol=imgReader.nrOfCol();
178 int ncropcol=0;
179 int ncroprow=0;
180 double dx=dx_opt[0];
181 double dy=dy_opt[0];
182 if(!dx||!dy){
183 dx=imgReader.getDeltaX();
184 dy=imgReader.getDeltaY();
185 }
186 //bounding box of cropped image
187 double cropulx=ulx_opt[0];
188 double cropuly=uly_opt[0];
189 double croplrx=lrx_opt[0];
190 double croplry=lry_opt[0];
191 //get bounding box from extentReader if defined
192 ImgReaderOgr extentReader;
193 if(extent_opt[0]!=""){
194 for(int iextent=0;iextent<extent_opt.size();++iextent){
195 extentReader.open(extent_opt[iextent]);
196 if(!(extentReader.getExtent(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
197 cerr << "Error: could not get extent from " << extent_opt[0] << endl;
198 exit(1);
199 }
200 if(!iextent){
201 cropulx=ulx_opt[0];
202 cropuly=uly_opt[0];
203 croplry=lry_opt[0];
204 croplrx=lrx_opt[0];
205 }
206 else{
207 if(ulx_opt[0]<cropulx)
208 cropulx=ulx_opt[0];
209 if(uly_opt[0]>cropuly)
210 cropuly=uly_opt[0];
211 if(lry_opt[0]<croplry)
212 croplry=lry_opt[0];
213 if(lrx_opt[0]>croplrx)
214 croplrx=lrx_opt[0];
215 }
216 extentReader.close();
217 }
218 }
219 double uli,ulj,lri,lrj;//image coordinates
220 double magicX=1,magicY=1;
221 if(ulx_opt[0]>=lrx_opt[0]){//default bounding box: no cropping
222 uli=0;
223 lri=imgReader.nrOfCol()-1;
224 ulj=0;
225 lrj=imgReader.nrOfRow()-1;
226 ncropcol=imgReader.nrOfCol();
227 ncroprow=imgReader.nrOfRow();
228 imgReader.getBoundingBox(cropulx,cropuly,croplrx,croplry);
229 imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
230 imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
231 ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx)));
232 ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy)));
233 }
234 else{
235 imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
236 imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
237
238 ncropcol=abs(static_cast<int>(ceil((croplrx-cropulx)/dx)));
239 ncroprow=abs(static_cast<int>(ceil((cropuly-croplry)/dy)));
240 uli=floor(uli);
241 ulj=floor(ulj);
242 lri=floor(lri);
243 lrj=floor(lrj);
244 }
245 double startCol=uli;
246 double endCol=lri;
247 double dcropcol=0;
248 double dcroprow=0;
249 double deltaX=imgReader.getDeltaX();
250 double deltaY=imgReader.getDeltaY();
251 dcropcol=(lri-uli+1)/(dx/deltaX);
252 dcroprow=(lrj-ulj+1)/(dy/deltaY);
253 if(verbose_opt[0]){
254 cout << "cropulx: " << cropulx << endl;
255 cout << "cropuly: " << cropuly << endl;
256 cout << "croplrx: " << croplrx << endl;
257 cout << "croplry: " << croplry << endl;
258 cout << "ncropcol: " << ncropcol << endl;
259 cout << "ncroprow: " << ncroprow << endl;
260 cout << "cropulx+ncropcol*dx: " << cropulx+ncropcol*dx << endl;
261 cout << "cropuly-ncroprow*dy: " << cropuly-ncroprow*dy << endl;
262 cout << "selected upper left column in input image: " << uli << endl;
263 cout << "selected upper left row of input image: " << ulj << endl;
264 cout << "selected lower right column in input image: " << lri << endl;
265 cout << "selected lower right row in input image: " << lrj << endl;
266 }
267 double gt[6];
268 imgReader.getGeoTransform(gt);
269 // imgWriter.setGeoTransform(gt);
270 // imgWriter.setProjection(imgReader.getProjection());
271
272 double readRow=0;
273 double readCol=0;
274 double lowerCol=0;
275 double upperCol=0;
276 int readncol=endCol-startCol+1;
277 //test
278 // vector<double> readBuffer(readncol);
279 vector<double> readBuffer(readncol+1);
280 // assert(imgWriter.isGeoRef());
281 if(band_opt.empty()){
282 for(int iband=0;iband<imgReader.nrOfBand();++iband)
283 band_opt.push_back(iband);
284 }
285 for(int iband=0;iband<band_opt.size();++iband){
286 assert(band_opt[iband]>=0);
287 assert(band_opt[iband]<imgReader.nrOfBand());
288 for(int irow=0;irow<ncroprow;++irow){
289 if(verbose_opt[0])
290 std::cout << irow << std::endl;
291 double x=0;
292 double y=0;
293 //convert irow to geo
294 // imgWriter.image2geo(0,irow,x,y);
295 imgReader.image2geo(0,irow,x,y);
296 //lookup corresponding row for irow in this file
297 imgReader.geo2image(x,y,readCol,readRow);
298 if(readRow<0||readRow>=imgReader.nrOfRow()){
299 if(verbose_opt[0]>1)
300 std::cout << "skipping row " << readRow << std::endl;
301 continue;
302 }
303 try{
304 //test
305 // imgReader.readData(readBuffer,startCol,endCol,readRow,band_opt[0],theResample);
306 if(endCol<imgReader.nrOfCol()-1)
307 imgReader.readData(readBuffer,startCol,endCol+1,readRow,band_opt[iband],theResample);
308 else
309 imgReader.readData(readBuffer,startCol,endCol,readRow,band_opt[iband],theResample);
310 for(int ib=0;ib<ncropcol;++ib){
311 // assert(imgWriter.image2geo(ib,irow,x,y));
312 assert(imgReader.image2geo(ib,irow,x,y));
313 //lookup corresponding row for irow in this file
314 imgReader.geo2image(x,y,readCol,readRow);
315 if(readCol<0||readCol>=imgReader.nrOfCol()){
316 if(oformat_opt[0]=="matrix"){
317 if(output_opt[0].empty())
318 std::cout << dstnodata_opt[0] << " ";
319 else
320 outputStream << dstnodata_opt[0] << " ";
321 }
322 else{
323 if(output_opt[0].empty())
324 std::cout << x << " " << y << " " << dstnodata_opt[0] << endl;
325 else
326 outputStream << x << " " << y << " " << dstnodata_opt[0] << endl;
327 }
328 }
329 else{
330 switch(theResample){
331 case(BILINEAR):
332 lowerCol=readCol-0.5;
333 lowerCol=static_cast<int>(lowerCol);
334 upperCol=readCol+0.5;
335 upperCol=static_cast<int>(upperCol);
336 if(lowerCol<0)
337 lowerCol=0;
338 if(upperCol>=imgReader.nrOfCol())
339 upperCol=imgReader.nrOfCol()-1;
340 if(oformat_opt[0]=="matrix"){
341 if(output_opt[0].empty())
342 std::cout << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " ";
343 else
344 outputStream << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " ";
345 }
346 else if(!imgReader.isNoData(readBuffer[upperCol-startCol])&&!imgReader.isNoData(readBuffer[lowerCol-startCol])){
347 if(output_opt[0].empty())
348 std::cout << x << " " << y << " " << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " " << endl;
349 else
350 outputStream << x << " " << y << " " << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " " << endl;
351 }
352 break;
353 default:
354 readCol=static_cast<int>(readCol);
355 readCol-=startCol;//we only start reading from startCol
356 assert(readCol>=0);
357 if(readCol>=readBuffer.size())
358 std::cout << "Error: " << readCol << " >= " << readBuffer.size() << std::endl;
359 assert(readCol<readBuffer.size());
360 if(oformat_opt[0]=="matrix"){
361 if(output_opt[0].empty())
362 std::cout << readBuffer[readCol] << " ";
363 else
364 outputStream << readBuffer[readCol] << " ";
365 }
366 else if(!imgReader.isNoData(readBuffer[readCol])){
367 if(output_opt[0].empty())
368 std::cout << x << " " << y << " " << readBuffer[readCol] << std::endl;
369 else
370 outputStream << x << " " << y << " " << readBuffer[readCol] << std::endl;
371 }
372 break;
373 }
374 }
375 }
376 }
377 catch(string errorstring){
378 cout << errorstring << endl;
379 exit(1);
380 }
381 if(oformat_opt[0]=="matrix"){
382 if(output_opt[0].empty())
383 std::cout << std::endl;
384 else
385 outputStream << std::endl;
386 }
387 }
388 }
389 if(extent_opt[0]!=""){
390 extentReader.close();
391 }
392 imgReader.close();
393 if(!output_opt[0].empty())
394 outputStream.close();
395}
ImgReaderGdal
Definition: ImgReaderGdal.h:37
ImgReaderOgr
Definition: ImgReaderOgr.h:36
Optionpk
Definition: Optionpk.h:107