pktools 2.6.7
Processing Kernel for geospatial data
pkinfo.cc
1/**********************************************************************
2pkinfo.cc: Report basic information from raster datasets (similar to gdalinfo)
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along 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/******************************************************************************/
90using namespace std;
91
92int 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}
Definition: Egcs.h:27
double getDeltaY(void) const
Get the pixel cell spacing in y.
int nrOfRow(void) const
Get the number of rows of this dataset.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
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)
bool getCenterPos(double &x, double &y) const
Get the center position of this dataset in georeferenced coordinates.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
char ** getMetadata()
Get the metadata of this dataset.
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
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.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
double getDeltaX(void) const
Get the pixel cell spacing in x.
std::string getInterleave() const
Get the band coding (interleave)
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
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...
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked,...
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
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
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...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
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...