| 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}
Egcs
Definition: Egcs.h:27
ImgRasterGdal::getDeltaY
double getDeltaY(void) const
Get the pixel cell spacing in y.
Definition: ImgRasterGdal.h:154
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::isGeoRef
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
Definition: ImgRasterGdal.h:104
ImgRasterGdal::getDriverDescription
std::string getDriverDescription() const
Get the GDAL driver description of this dataset.
Definition: ImgRasterGdal.cc:131
ImgRasterGdal::geo2image
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
Definition: ImgRasterGdal.cc:387
ImgRasterGdal::getCenterPos
bool getCenterPos(double &x, double &y) const
Get the center position of this dataset in georeferenced coordinates.
Definition: ImgRasterGdal.cc:359
ImgRasterGdal::getProjection
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
Definition: ImgRasterGdal.cc:40
ImgRasterGdal::getMetadata
char ** getMetadata()
Get the metadata of this dataset.
Definition: ImgRasterGdal.cc:214
ImgRasterGdal::getGeoTransform
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
Definition: ImgRasterGdal.cc:200
ImgRasterGdal::getColorTable
GDALColorTable * getColorTable(int band=0) const
Get the GDAL color table for this dataset as an instance of the GDALColorTable class.
Definition: ImgRasterGdal.cc:118
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgRasterGdal::getDeltaX
double getDeltaX(void) const
Get the pixel cell spacing in x.
Definition: ImgRasterGdal.h:152
ImgRasterGdal::getInterleave
std::string getInterleave() const
Get the band coding (interleave)
Definition: ImgRasterGdal.cc:293
ImgRasterGdal::getRasterBand
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
Definition: ImgRasterGdal.cc:105
ImgRasterGdal::pushNoDataValue
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
Definition: ImgRasterGdal.cc:492
ImgRasterGdal::GDALSetNoDataValue
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...
Definition: ImgRasterGdal.h:142
ImgRasterGdal::getDataType
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
Definition: ImgRasterGdal.cc:92
ImgRasterGdal::covers
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked,...
Definition: ImgRasterGdal.cc:443
ImgRasterGdal::getBoundingBox
bool getBoundingBox(double &ulx, double &uly, double &lrx, double &lry) const
Get the bounding box of this dataset in georeferenced coordinates.
Definition: ImgRasterGdal.cc:333
ImgReaderGdal
Definition: ImgReaderGdal.h:37
ImgReaderGdal::readData
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
ImgReaderGdal::getMinMax
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...
Definition: ImgReaderGdal.cc:161
ImgReaderGdal::close
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
Definition: ImgReaderGdal.cc:42
ImgReaderGdal::open
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
Definition: ImgReaderGdal.cc:36
ImgReaderGdal::getRefPix
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...
Definition: ImgReaderGdal.cc:398
ImgReaderOgr
Definition: ImgReaderOgr.h:36
Optionpk
Definition: Optionpk.h:107