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 }
STL namespace.