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 }
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.
STL namespace.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
std::string getDriverDescription() const
Get the GDAL driver description of this dataset.
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
int nrOfRow(void) const
Get the number of rows of this dataset.
bool getCenterPos(double &x, double &y) const
Get the center position of this dataset in georeferenced coordinates.
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
double getDeltaX(void) const
Get the pixel cell spacing in x.
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
char ** getMetadata()
Get the metadata of this dataset.
double getDeltaY(void) const
Get the pixel cell spacing in y.
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...
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
Definition: Egcs.h:26
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...
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
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.
int nrOfBand(void) const
Get the number of bands of this dataset.
std::string getInterleave() const
Get the band coding (interleave)