pktools  2.6.7
Processing Kernel for geospatial data
pkinfo.cc
1 /**********************************************************************
2 pkinfo.cc: Report basic information from raster datasets (similar to gdalinfo)
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 <sstream>
21 #include <list>
22 #include "base/Optionpk.h"
23 #include "algorithms/Egcs.h"
24 #include "imageclasses/ImgReaderGdal.h"
25 #include "imageclasses/ImgReaderOgr.h"
26 
27 /******************************************************************************/
90 using namespace std;
91 
92 int main(int argc, char *argv[])
93 {
94  Optionpk<std::string> input_opt("i","input","Input image file");
95  Optionpk<bool> bbox_opt("bb", "bbox", "Shows bounding box ", false,0);
96  Optionpk<bool> bbox_te_opt("te", "te", "Shows bounding box in GDAL format: xmin ymin xmax ymax ", false,0);
97  Optionpk<bool> center_opt("c", "center", "Image center in projected X,Y coordinates ", false,0);
98  Optionpk<bool> colorTable_opt("ct", "colortable", "Shows colour table ", false,0);
99  Optionpk<bool> samples_opt("ns", "nsample", "Number of samples in image ", false,0);
100  Optionpk<bool> lines_opt("nl", "nline", "Number of lines in image ", false,0);
101  Optionpk<bool> nband_opt("nb", "nband", "Show number of bands in image", false,0);
102  Optionpk<short> band_opt("b", "band", "Band specific information", 0,0);
103  Optionpk<bool> dx_opt("dx", "dx", "Gets resolution in x (in m)", false,0);
104  Optionpk<bool> dy_opt("dy", "dy", "Gets resolution in y (in m)", false,0);
105  Optionpk<bool> minmax_opt("mm", "minmax", "Shows min and max value of the image ", false,0);
106  Optionpk<bool> min_opt("min", "minimum", "Shows min value of the image ", false,0);
107  Optionpk<bool> max_opt("max", "maximum", "Shows max value of the image ", false,0);
108  Optionpk<bool> stat_opt("stats", "statistics", "Shows statistics (min,max, mean and stdDev of the image)", false,0);
109  Optionpk<bool> projection_opt("a_srs", "a_srs", "Shows projection of the image ", false,0);
110  Optionpk<bool> geo_opt("geo", "geo", "Gets geotransform ", false,0);
111  Optionpk<bool> interleave_opt("il", "interleave", "Shows interleave ", false,0);
112  Optionpk<bool> filename_opt("f", "filename", "Shows image filename ", false,0);
113  Optionpk<bool> cover_opt("cover", "cover", "Print filename to stdout if current image covers the provided coordinates via bounding box, (x y) coordinates or extent of vector file", false,0);
114  Optionpk<double> x_opt("x", "xpos", "x pos");
115  Optionpk<double> y_opt("y", "ypos", "y pos");
116  Optionpk<bool> read_opt("r", "read", "Reads row y (in projected coordinates if geo option is set, otherwise in image coordinates, 0 based)",false,0);
117  Optionpk<bool> refpixel_opt("ref", "reference", "Gets reference pixel (lower left corner of center of gravity pixel)", false,0);
118  Optionpk<bool> driver_opt("of", "oformat", "Gets driver description ", false,0);
119  Optionpk<std::string> extent_opt("e", "extent", "Gets boundary from vector file");
120  Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box");
121  Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box");
122  Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box");
123  Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box");
124  Optionpk<bool> type_opt("ot", "otype", "Returns data type", false,0);
125  Optionpk<bool> description_opt("d", "description", "Returns image description", false,0);
126  Optionpk<bool> metadata_opt("meta", "meta", "Shows meta data ", false,0);
127  Optionpk<double> nodata_opt("nodata", "nodata", "Sets no data value(s) for calculations (nodata values in input image)");
128 
129  bool doProcess;//stop process when program was invoked with help option (-h --help)
130  try{
131  doProcess=input_opt.retrieveOption(argc,argv);
132  bbox_opt.retrieveOption(argc,argv);
133  bbox_te_opt.retrieveOption(argc,argv);
134  center_opt.retrieveOption(argc,argv);
135  colorTable_opt.retrieveOption(argc,argv);
136  samples_opt.retrieveOption(argc,argv);
137  lines_opt.retrieveOption(argc,argv);
138  nband_opt.retrieveOption(argc,argv);
139  band_opt.retrieveOption(argc,argv);
140  dx_opt.retrieveOption(argc,argv);
141  dy_opt.retrieveOption(argc,argv);
142  minmax_opt.retrieveOption(argc,argv);
143  min_opt.retrieveOption(argc,argv);
144  max_opt.retrieveOption(argc,argv);
145  stat_opt.retrieveOption(argc,argv);
146  projection_opt.retrieveOption(argc,argv);
147  geo_opt.retrieveOption(argc,argv);
148  interleave_opt.retrieveOption(argc,argv);
149  filename_opt.retrieveOption(argc,argv);
150  cover_opt.retrieveOption(argc,argv);
151  x_opt.retrieveOption(argc,argv);
152  y_opt.retrieveOption(argc,argv);
153  read_opt.retrieveOption(argc,argv);
154  refpixel_opt.retrieveOption(argc,argv);
155  driver_opt.retrieveOption(argc,argv);
156  extent_opt.retrieveOption(argc,argv);
157  ulx_opt.retrieveOption(argc,argv);
158  uly_opt.retrieveOption(argc,argv);
159  lrx_opt.retrieveOption(argc,argv);
160  lry_opt.retrieveOption(argc,argv);
161  type_opt.retrieveOption(argc,argv);
162  description_opt.retrieveOption(argc,argv);
163  metadata_opt.retrieveOption(argc,argv);
164  nodata_opt.retrieveOption(argc,argv);
165  }
166  catch(std::string predefinedString){
167  std::cout << predefinedString << std::endl;
168  exit(0);
169  }
170  if(!doProcess){
171  cout << endl;
172  cout << "Usage: pkinfo -i input [options]" << endl;
173  cout << endl;
174  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
175  exit(0);//help was invoked, stop processing
176  }
177  //for union
178  double maxLRX=0;
179  double maxULY=0;
180  double minULX=0;
181  double minLRY=0;
182 
183  //for intersect
184  double minLRX=0;
185  double minULY=0;
186  double maxULX=0;
187  double maxLRY=0;
188 
189  double theULX, theULY, theLRX, theLRY;
190  //get bounding box from extentReader if defined
191  ImgReaderOgr extentReader;
192  if(extent_opt.size()){
193  extentReader.open(extent_opt[0]);
194  if(!(extentReader.getExtent(theULX,theULY, theLRX, theLRY))){
195  std::cerr << "Error: could not get extent from " << extent_opt[0] << std::endl;
196  exit(1);
197  }
198  ulx_opt.push_back(theULX);
199  uly_opt.push_back(theULY);
200  lrx_opt.push_back(theLRX);
201  lry_opt.push_back(theLRY);
202  if(input_opt.empty()){//report bounding box from extent file instead
203  if(bbox_te_opt[0])
204  std::cout << std::setprecision(12) << "-te " << theULX << " " << theLRY << " " << theLRX << " " << theULY;
205  else
206  std::cout << std::setprecision(12) << "--ulx=" << theULX << " --uly=" << theULY << " --lrx=" << theLRX << " --lry=" << theLRY << " ";
207  }
208  }
209 
210  ImgReaderGdal imgReader;
211  for(int ifile=0;ifile<input_opt.size();++ifile){
212  imgReader.open(input_opt[ifile]);
213  for(int inodata=0;inodata<nodata_opt.size();++inodata){
214  if(!inodata)
215  imgReader.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
216  imgReader.pushNoDataValue(nodata_opt[inodata]);
217  }
218  if(filename_opt[0])
219  std::cout << " --input " << input_opt[ifile] << " ";
220  if(center_opt[0]){
221  double theX, theY;
222  imgReader.getCenterPos(theX,theY);
223  std::cout << std::setprecision(12) << " -x " << theX << " -y " << theY << " ";
224  }
225  if(refpixel_opt[0]){
226  assert(band_opt[0]<imgReader.nrOfBand());
227  Egcs egcs;
228  double refX,refY;
229  //get center of reference (center of gravity) pixel in image
230  imgReader.getRefPix(refX,refY,band_opt[0]);
231  std::cout << std::setprecision(12) << "-x " << refX << " -y " << refY << std::endl;
232  egcs.setLevel(egcs.res2level(imgReader.getDeltaX()));
233  // unsigned short theLevel=egcs.getLevel(imgReader.getDeltaX());
234  // egcs.setLevel(theLevel);
235  //cout << "cell code at level " << egcs.getLevel() << " (resolution is " << egcs.getResolution() << "): " << egcs.geo2cell(refX,refY) << endl;
236  }
237  if(bbox_opt[0]||bbox_te_opt[0]){
238  imgReader.getBoundingBox(theULX,theULY,theLRX,theLRY);
239  if(bbox_te_opt[0])
240  std::cout << std::setprecision(12) << "-te " << theULX << " " << theLRY << " " << theLRX << " " << theULY;
241  else
242  std::cout << std::setprecision(12) << "--ulx=" << theULX << " --uly=" << theULY << " --lrx=" << theLRX << " --lry=" << theLRY << " ";
243  if(!ifile){
244  maxLRX=theLRX;
245  maxULY=theULY;
246  minULX=theULX;
247  minLRY=theLRY;
248 
249  minLRX=theLRX;
250  minULY=theULY;
251  maxULX=theULX;
252  maxLRY=theLRY;
253  }
254  else{
255  maxLRX=(theLRX>maxLRX)?theLRX:maxLRX;
256  maxULY=(theULY>maxULY)?theULY:maxULY;
257  minULX=(theULX<minULX)?theULX:minULX;
258  minLRY=(theLRY<minLRY)?theLRY:minLRY;
259 
260  minLRX=(theLRX<minLRX)?theLRX:minLRX;
261  minULY=(theULY<minULY)?theULY:minULY;
262  maxULX=(theULX>maxULX)?theULX:maxULX;
263  maxLRY=(theLRY>maxLRY)?theLRY:maxLRY;
264  }
265  }
266  if(dx_opt[0])
267  std::cout << "--dx " << imgReader.getDeltaX() << " ";
268  if(dy_opt[0])
269  std::cout << "--dy " << imgReader.getDeltaY() << " ";
270  if(cover_opt[0]){
271  if(ulx_opt.size()&&uly_opt.size()&&lrx_opt.size()&&lry_opt.size()){
272  if(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))
273  std::cout << " -i " << input_opt[ifile] << " ";
274  }
275  else if(x_opt.size()&&y_opt.size()){
276  if(imgReader.covers(x_opt[0],y_opt[0]))
277  std::cout << " -i " << input_opt[ifile] << " ";
278  }
279  else{
280  std::cerr << "Error: failing extent (-e), bounding box or x and y position to define coverage" << std::endl;
281  exit(1);
282  }
283  }
284  else if(ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size()){
285  double ulx,uly,lrx,lry;
286  imgReader.getBoundingBox(ulx,uly,lrx,lry);
287  if(ulx_opt.size())
288  std::cout << " --ulx=" << std::fixed << ulx << " ";
289  if(uly_opt.size())
290  std::cout << " --uly=" << std::fixed << uly << " ";
291  if(lrx_opt.size())
292  std::cout << " --lrx=" << std::fixed << lrx << " ";
293  if(lry_opt.size())
294  std::cout << " --lry=" << std::fixed << lry << " ";
295  }
296  if(colorTable_opt[0]){
297  GDALColorTable* colorTable=imgReader.getColorTable();
298  if(colorTable!=NULL){
299  for(int index=0;index<colorTable->GetColorEntryCount();++index){
300  GDALColorEntry sEntry=*(colorTable->GetColorEntry(index));
301  std::cout << index << " " << sEntry.c1 << " " << sEntry.c2 << " " << sEntry.c3 << " " << sEntry.c4 << std::endl;
302  }
303  }
304  else
305  std::cout << "-ct none ";
306  }
307  if(samples_opt[0])
308  std::cout << "--nsample " << imgReader.nrOfCol() << " ";
309  if(lines_opt[0])
310  std::cout << "--nline " << imgReader.nrOfRow() << " ";
311  if(nband_opt[0])
312  std::cout << "--nband " << imgReader.nrOfBand() << " ";
313  double minValue=0;
314  double maxValue=0;
315  double meanValue=0;
316  double stdDev=0;
317  int nband=band_opt.size();
318  if(band_opt[0]<0)
319  nband=imgReader.nrOfBand();
320  for(int iband=0;iband<nband;++iband){
321  unsigned short theBand=(band_opt[0]<0)? iband : band_opt[iband];
322  if(stat_opt[0]){
323  assert(theBand<imgReader.nrOfBand());
324  GDALProgressFunc pfnProgress;
325  void* pProgressData;
326  GDALRasterBand* rasterBand;
327  rasterBand=imgReader.getRasterBand(theBand);
328  rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
329  std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " ";
330  }
331 
332  if(minmax_opt[0]||min_opt[0]||max_opt[0]){
333  assert(theBand<imgReader.nrOfBand());
334  if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
335  double uli,ulj,lri,lrj;
336  imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
337  imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
338  imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),theBand,minValue,maxValue);
339  }
340  else
341  imgReader.getMinMax(minValue,maxValue,theBand);
342  if(minmax_opt[0])
343  std::cout << "-min " << minValue << " -max " << maxValue << " ";
344  else{
345  if(min_opt[0])
346  std::cout << "-min " << minValue << " ";
347  if(max_opt[0])
348  std::cout << "-max " << maxValue << " ";
349  }
350  }
351  }
352  if(projection_opt[0]){
353  if(imgReader.isGeoRef())
354  std::cout << " -a_srs " << imgReader.getProjection() << " ";
355  else
356  std::cout << " -a_srs none" << " ";
357  }
358  if(geo_opt[0]&&!read_opt[0]){
359  std::cout << " -geo " << std::setprecision(12) << imgReader.getGeoTransform();
360  }
361  if(interleave_opt[0]){
362  std::cout << " --interleave " << imgReader.getInterleave() << " ";
363  }
364  if(type_opt[0]){
365  std::cout << "--otype " << GDALGetDataTypeName(imgReader.getDataType(band_opt[0])) << " ";
366  // std::cout << " -ot " << GDALGetDataTypeName(imgReader.getDataType(band_opt[0])) << " (" << static_cast<short>(imgReader.getDataType(band_opt[0])) << ")" << std::endl;
367  }
368  if(description_opt[0]){
369  // try{
370  // std::cout << "image description: " << imgReader.getImageDescription() << std::endl;
371  // }
372  // catch(...){
373  // std::cout << "caught" << std::endl;
374  // }
375  std::list<std::string> metaData;
376  imgReader.getMetadata(metaData);
377  std::list<std::string>::const_iterator lit=metaData.begin();
378  std::cout << " --description ";
379  while(lit!=metaData.end())
380  std::cout << *(lit++) << " ";
381  }
382  if(metadata_opt[0]){
383  std::cout << "Metadata: " << std::endl;
384  std::list<std::string> lmeta;
385  imgReader.getMetadata(lmeta);
386  std::list<std::string>::const_iterator lit=lmeta.begin();
387  while(lit!=lmeta.end()){
388  std::cout << *lit << std::endl;
389  ++lit;
390  }
391 // char** cmetadata=imgReader.getMetadata();
392 // while(*cmetadata!=NULL){
393 // std::cout << *(cmetadata) << std::endl;
394 // ++cmetadata;
395 // }
396  }
397  if(read_opt[0]){
398  // int nband=band_opt.size();
399  // if(band_opt[0]<0)
400  // nband=imgReader.nrOfBand();
401  std::cout.precision(12);
402  for(int iband=0;iband<nband;++iband){
403  unsigned short theBand=(band_opt[0]<0)? iband : band_opt[iband];
404  std::vector<float> rowBuffer;//buffer will be resized in readdata
405  for(int iy=0;iy<y_opt.size();++iy){
406  double theRow=y_opt[iy];
407  int ncol=(x_opt.size())? x_opt.size() : imgReader.nrOfCol();
408  for(int ix=0;ix<ncol;++ix){
409  double theCol=ix;
410  if(x_opt.size()){
411  if(geo_opt[0])
412  imgReader.geo2image(x_opt[ix],y_opt[iy],theCol,theRow);
413  else
414  theCol=x_opt[ix];
415  }
416  assert(theRow>=0);
417  assert(theRow<imgReader.nrOfRow());
418  imgReader.readData(rowBuffer, static_cast<int>(theRow), theBand);
419  assert(theCol<rowBuffer.size());
420  std::cout << rowBuffer[static_cast<int>(theCol)] << " ";
421  }
422  std::cout << std::endl;
423  }
424  }
425  }
426  if(driver_opt[0])
427  std::cout << " --oformat " << imgReader.getDriverDescription() << " ";
428  imgReader.close();
429  }
430  if((bbox_opt[0]||bbox_te_opt[0])&&input_opt.size()>1){
431  if(bbox_te_opt[0])
432  std::cout << std::setprecision(12) << "-te " << minULX << " " << minLRY << " " << maxLRX << " " << maxULY;
433  else
434  std::cout << "union bounding box: " << std::setprecision(12) << "--ulx=" << minULX << " --uly=" << maxULY << " --lrx=" << maxLRX << " --lry=" << minLRY << std::endl;
435  if(maxULX<minLRX&&minULY>maxLRY){
436  if(bbox_te_opt[0])
437  std::cout << "intersect bounding box: " << std::setprecision(12) << "-te " << maxULX << " " << maxLRY << " " << minLRX << " --lry=" << minULY << std::endl;
438  else
439  std::cout << "intersect bounding box: " << std::setprecision(12) << "--ulx=" << maxULX << " --uly=" << minULY << " --lrx=" << minLRX << " --lry=" << maxLRY << std::endl;
440  }
441  else
442  std::cout << "no intersect" << std::endl;
443  }
444  if(!read_opt[0])
445  std::cout << std::endl;
446 }
std::string getDriverDescription() const
Get the GDAL driver description of this dataset.
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.
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) ?
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.
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked, irrespective of no data values.
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.
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
char ** getMetadata()
Get the metadata of this dataset.
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 getDeltaX(void) const
Get the pixel cell spacing in x.
bool getCenterPos(double &x, double &y) const
Get the center position of this dataset in georeferenced coordinates.
std::string getInterleave() const
Get the band coding (interleave)
Definition: Egcs.h:26
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
void getRefPix(double &refX, double &refY, int band=0)
Calculate the reference pixel as the centre of gravity pixel (weighted average of all values not taki...
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
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.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98