| pktools 2.6.7 Processing Kernel for geospatial data

ImgReaderGdal.cc
1/**********************************************************************
2 ImgReaderGdal.cc: class to read raster files using GDAL API library
3 Copyright (C) 2008-2016 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 <assert.h>
21#include <sstream>
22#include <iostream>
23#include <gsl/gsl_cdf.h>
24#include "cpl_conv.h" // for CPLMalloc()
25#include "ImgReaderGdal.h"
26
27 ImgReaderGdal::ImgReaderGdal(void){};
28
29 ImgReaderGdal::~ImgReaderGdal(void){};
30
36 void ImgReaderGdal::open(const std::string& filename, const GDALAccess& readMode)
37{
38 m_filename = filename;
39 setCodec(readMode);
40}
41
42 void ImgReaderGdal::close(void)
43{
44 ImgRasterGdal::close();
45}
46
50 void ImgReaderGdal::setCodec(const GDALAccess& readMode)
51{
52 GDALAllRegister();
53 // m_gds = (GDALDataset *) GDALOpen(m_filename.c_str(), readMode );
54#if GDAL_VERSION_MAJOR < 2
55 GDALAllRegister();
56 m_gds = (GDALDataset *) GDALOpen(m_filename.c_str(), readMode );
57#else
58 GDALAllRegister();
59 if(readMode==GA_ReadOnly)
60 m_gds = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_READONLY|GDAL_OF_RASTER, NULL, NULL, NULL);
61 else if(readMode==GA_Update)
62 m_gds = (GDALDataset*) GDALOpenEx(m_filename.c_str(), GDAL_OF_UPDATE|GDAL_OF_RASTER, NULL, NULL, NULL);
63#endif
64
65 if(m_gds == NULL){
66 std::string errorString="FileOpenError";
67 throw(errorString);
68 }
69 m_ncol= m_gds->GetRasterXSize();
70 m_nrow= m_gds->GetRasterYSize();
71 m_nband= m_gds->GetRasterCount();
72 double adfGeoTransform[6];
73 m_gds->GetGeoTransform( adfGeoTransform );
74 m_gt[0]=adfGeoTransform[0];
75 m_gt[1]=adfGeoTransform[1];
76 m_gt[2]=adfGeoTransform[2];
77 m_gt[3]=adfGeoTransform[3];
78 m_gt[4]=adfGeoTransform[4];
79 m_gt[5]=adfGeoTransform[5];
80 m_projection=m_gds->GetProjectionRef();
81}
82
89 double ImgReaderGdal::getMin(int& x, int& y, int band){
90 double minValue=0;
91 std::vector<double> lineBuffer(nrOfCol());
92 bool isValid=false;
93 for(int irow=0;irow<nrOfRow();++irow){
94 readData(lineBuffer,irow,band);
95 for(int icol=0;icol<nrOfCol();++icol){
96 if(isNoData(lineBuffer[icol]))
97 continue;
98 if(isValid){
99 if(lineBuffer[icol]<minValue){
100 y=irow;
101 x=icol;
102 minValue=lineBuffer[icol];
103 }
104 }
105 else{
106 y=irow;
107 x=icol;
108 minValue=lineBuffer[icol];
109 isValid=true;
110 }
111 }
112 }
113 if(isValid)
114 return minValue;
115 else
116 throw(static_cast<std::string>("Warning: not initialized"));
117}
118
125 double ImgReaderGdal::getMax(int& x, int& y, int band){
126 double maxValue=0;
127 std::vector<double> lineBuffer(nrOfCol());
128 bool isValid=false;
129 for(int irow=0;irow<nrOfRow();++irow){
130 readData(lineBuffer,irow,band);
131 for(int icol=0;icol<nrOfCol();++icol){
132 if(isNoData(lineBuffer[icol]))
133 continue;
134 if(isValid){
135 if(lineBuffer[icol]>maxValue){
136 y=irow;
137 x=icol;
138 maxValue=lineBuffer[icol];
139 }
140 }
141 else{
142 y=irow;
143 x=icol;
144 maxValue=lineBuffer[icol];
145 isValid=true;
146 }
147 }
148 }
149 if(isValid)
150 return maxValue;
151 else
152 throw(static_cast<std::string>("Warning: not initialized"));
153}
154
161 void ImgReaderGdal::getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double& minValue, double& maxValue)
162{
163 bool isConstraint=(maxValue>minValue);
164 double minConstraint=minValue;
165 double maxConstraint=maxValue;
166 std::vector<double> lineBuffer(endCol-startCol+1);
167 bool isValid=false;
168 assert(endRow<nrOfRow());
169 for(int irow=startCol;irow<endRow+1;++irow){
170 readData(lineBuffer,startCol,endCol,irow,band);
171 for(int icol=0;icol<lineBuffer.size();++icol){
172 if(isNoData(lineBuffer[icol]))
173 continue;
174 if(isValid){
175 if(isConstraint){
176 if(lineBuffer[icol]<minConstraint)
177 continue;
178 if(lineBuffer[icol]>maxConstraint)
179 continue;
180 }
181 if(lineBuffer[icol]<minValue)
182 minValue=lineBuffer[icol];
183 if(lineBuffer[icol]>maxValue)
184 maxValue=lineBuffer[icol];
185 }
186 else{
187 if(isConstraint){
188 if(lineBuffer[icol]<minConstraint)
189 continue;
190 if(lineBuffer[icol]>maxConstraint)
191 continue;
192 }
193 minValue=lineBuffer[icol];
194 maxValue=lineBuffer[icol];
195 isValid=true;
196 }
197 }
198 }
199 if(!isValid)
200 throw(static_cast<std::string>("Warning: not initialized"));
201}
202
208 void ImgReaderGdal::getMinMax(double& minValue, double& maxValue, int band)
209{
210 bool isConstraint=(maxValue>minValue);
211 double minConstraint=minValue;
212 double maxConstraint=maxValue;
213 std::vector<double> lineBuffer(nrOfCol());
214 bool isValid=false;
215 for(int irow=0;irow<nrOfRow();++irow){
216 readData(lineBuffer,irow,band);
217 for(int icol=0;icol<nrOfCol();++icol){
218 if(isNoData(lineBuffer[icol]))
219 continue;
220 if(isValid){
221 if(isConstraint){
222 if(lineBuffer[icol]<minConstraint)
223 continue;
224 if(lineBuffer[icol]>maxConstraint)
225 continue;
226 }
227 if(lineBuffer[icol]<minValue)
228 minValue=lineBuffer[icol];
229 if(lineBuffer[icol]>maxValue)
230 maxValue=lineBuffer[icol];
231 }
232 else{
233 if(isConstraint){
234 if(lineBuffer[icol]<minConstraint)
235 continue;
236 if(lineBuffer[icol]>maxConstraint)
237 continue;
238 }
239 minValue=lineBuffer[icol];
240 maxValue=lineBuffer[icol];
241 isValid=true;
242 }
243 }
244 }
245 if(!isValid)
246 throw(static_cast<std::string>("Warning: not initialized"));
247}
248
257 double ImgReaderGdal::getHistogram(std::vector<double>& histvector, double& min, double& max, unsigned int& nbin, int theBand, bool kde){
258 double minValue=0;
259 double maxValue=0;
260
261 if(min>=max)
262 getMinMax(minValue,maxValue,theBand);
263 else{
264 minValue=min;
265 maxValue=max;
266 }
267 if(min<max&&min>minValue)
268 minValue=min;
269 if(min<max&&max<maxValue)
270 maxValue=max;
271 min=minValue;
272 max=maxValue;
273
274 double sigma=0;
275 if(kde){
276 double meanValue=0;
277 double stdDev=0;
278 GDALProgressFunc pfnProgress;
279 void* pProgressData;
280 GDALRasterBand* rasterBand;
281 rasterBand=getRasterBand(theBand);
282 rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
283 //rest minvalue and MaxValue as ComputeStatistics does not account for nodata, scale and offset
284 minValue=min;
285 maxValue=max;
286
287 if(m_scale.size()>theBand){
288 stdDev*=m_scale[theBand];
289 }
290 sigma=1.06*stdDev*pow(getNvalid(theBand),-0.2);
291 }
292
293 double scale=0;
294 if(maxValue>minValue){
295 if(nbin==0)
296 nbin=maxValue-minValue+1;
297 scale=static_cast<double>(nbin-1)/(maxValue-minValue);
298 }
299 else
300 nbin=1;
301 assert(nbin>0);
302 if(histvector.size()!=nbin){
303 histvector.resize(nbin);
304 for(int i=0;i<nbin;histvector[i++]=0);
305 }
306 double nvalid=0;
307 unsigned long int nsample=0;
308 unsigned long int ninvalid=0;
309 std::vector<double> lineBuffer(nrOfCol());
310 for(int irow=0;irow<nrOfRow();++irow){
311 readData(lineBuffer,irow,theBand);
312 for(int icol=0;icol<nrOfCol();++icol){
313 if(isNoData(lineBuffer[icol]))
314 ++ninvalid;
315 else if(lineBuffer[icol]>maxValue)
316 ++ninvalid;
317 else if(lineBuffer[icol]<minValue)
318 ++ninvalid;
319 else if(nbin==1)
320 ++histvector[0];
321 else{//scale to [0:nbin]
322 if(sigma>0){
323 //create kde for Gaussian basis function
324 //todo: speed up by calculating first and last bin with non-zero contriubtion...
325 //todo: calculate real surface below pdf by using gsl_cdf_gaussian_P(x-mean+binsize,sigma)-gsl_cdf_gaussian_P(x-mean,sigma)
326 for(int ibin=0;ibin<nbin;++ibin){
327 double icenter=minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin;
328 double thePdf=gsl_ran_gaussian_pdf(lineBuffer[icol]-icenter, sigma);
329 histvector[ibin]+=thePdf;
330 nvalid+=thePdf;
331 }
332 }
333 else{
334 int theBin=static_cast<unsigned long int>(scale*(lineBuffer[icol]-minValue));
335 assert(theBin>=0);
336 assert(theBin<nbin);
337 ++histvector[theBin];
338 ++nvalid;
339 }
340 // else if(lineBuffer[icol]==maxValue)
341 // ++histvector[nbin-1];
342 // else
343 // ++histvector[static_cast<int>(static_cast<double>(lineBuffer[icol]-minValue)/(maxValue-minValue)*(nbin-1))];
344 }
345 }
346 }
347 // unsigned long int nvalid=nrOfCol()*nrOfRow()-ninvalid;
348 return nvalid;
349}
350
355 void ImgReaderGdal::getRange(std::vector<short>& range, int band)
356{
357 std::vector<short> lineBuffer(nrOfCol());
358 range.clear();
359 for(int irow=0;irow<nrOfRow();++irow){
360 readData(lineBuffer,irow,band);
361 for(int icol=0;icol<nrOfCol();++icol){
362 if(find(range.begin(),range.end(),lineBuffer[icol])==range.end())
363 range.push_back(lineBuffer[icol]);
364 }
365 }
366 sort(range.begin(),range.end());
367}
368
373 unsigned long int ImgReaderGdal::getNvalid(int band)
374{
375 unsigned long int nvalid=0;
376 if(m_noDataValues.size()){
377 std::vector<double> lineBuffer(nrOfCol());
378 for(int irow=0;irow<nrOfRow();++irow){
379 readData(lineBuffer,irow,band);
380 for(int icol=0;icol<nrOfCol();++icol){
381 if(isNoData(lineBuffer[icol]))
382 continue;
383 else
384 ++nvalid;
385 }
386 }
387 return nvalid;
388 }
389 else
390 return(nrOfCol()*nrOfRow());
391}
392
398 void ImgReaderGdal::getRefPix(double& refX, double &refY, int band)
399{
400 std::vector<double> lineBuffer(nrOfCol());
401 double validCol=0;
402 double validRow=0;
403 int nvalidCol=0;
404 int nvalidRow=0;
405 for(int irow=0;irow<nrOfRow();++irow){
406 readData(lineBuffer,irow,band);
407 for(int icol=0;icol<nrOfCol();++icol){
408 // bool valid=(find(m_noDataValues.begin(),m_noDataValues.end(),lineBuffer[icol])==m_noDataValues.end());
409 // if(valid){
410 if(!isNoData(lineBuffer[icol])){
411 validCol+=icol+1;
412 ++nvalidCol;
413 validRow+=irow+1;
414 ++nvalidRow;
415 }
416 }
417 }
418 if(isGeoRef()){
419 //reference coordinate is lower left corner of pixel in center of gravity
420 //we need geo coordinates for exactly this location: validCol(Row)/nvalidCol(Row)-0.5
421 double cgravi=validCol/nvalidCol-0.5;
422 double cgravj=validRow/nvalidRow-0.5;
423 double refpixeli=floor(cgravi);
424 double refpixelj=ceil(cgravj-1);
425 //but image2geo provides location at center of pixel (shifted half pixel right down)
426 image2geo(refpixeli,refpixelj,refX,refY);
427 //refX and refY now refer to center of gravity pixel
428 refX-=0.5*getDeltaX();//shift to left corner
429 refY-=0.5*getDeltaY();//shift to lower left corner
430 }
431 else{
432 refX=floor(validCol/nvalidCol-0.5);//left corner
433 refY=floor(validRow/nvalidRow-0.5);//upper corner
434 //shift to lower left corner of pixel
435 refY+=1;
436 }
437}
ImgRasterGdal::getDeltaY
double getDeltaY(void) const
Get the pixel cell spacing in y.
Definition: ImgRasterGdal.h:154
ImgRasterGdal::close
virtual void close(void)
Close the image.
Definition: ImgRasterGdal.cc:32
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::m_nband
int m_nband
number of bands in this dataset
Definition: ImgRasterGdal.h:211
ImgRasterGdal::m_ncol
int m_ncol
number of columns in this dataset
Definition: ImgRasterGdal.h:207
ImgRasterGdal::isNoData
bool isNoData(double value) const
Check if value is nodata in this dataset.
Definition: ImgRasterGdal.h:136
ImgRasterGdal::m_filename
std::string m_filename
filename of this dataset
Definition: ImgRasterGdal.h:203
ImgRasterGdal::m_gds
GDALDataset * m_gds
instance of the GDAL dataset of this dataset
Definition: ImgRasterGdal.h:205
ImgRasterGdal::m_nrow
int m_nrow
number of rows in this dataset
Definition: ImgRasterGdal.h:209
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::m_noDataValues
std::vector< double > m_noDataValues
no data values for this dataset
Definition: ImgRasterGdal.h:219
ImgRasterGdal::getRasterBand
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
Definition: ImgRasterGdal.cc:105
ImgRasterGdal::m_scale
std::vector< double > m_scale
Vector containing the scale factor to be applied (one scale value for each band)
Definition: ImgRasterGdal.h:221
ImgRasterGdal::image2geo
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
Definition: ImgRasterGdal.cc:416
ImgRasterGdal::m_gt
double m_gt[6]
geotransform information of this dataset
Definition: ImgRasterGdal.h:215
ImgReaderGdal::setCodec
void setCodec(const GDALAccess &readMode=GA_ReadOnly)
Set GDAL dataset number of columns, rows, bands and geotransform.
Definition: ImgReaderGdal.cc:50
ImgReaderGdal::getNvalid
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
Definition: ImgReaderGdal.cc:373
ImgReaderGdal::~ImgReaderGdal
~ImgReaderGdal(void)
destructor
Definition: ImgReaderGdal.cc:29
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::getMax
double getMax(int &col, int &row, int band=0)
Get the maximum cell values for a specific band and report the column and row in which the maximum va...
Definition: ImgReaderGdal.cc:125
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
ImgReaderGdal::ImgReaderGdal
ImgReaderGdal(void)
default constructor. Image needs to be opened later with one of the open methods.
Definition: ImgReaderGdal.cc:27
ImgReaderGdal::getRange
void getRange(std::vector< short > &range, int Band=0)
Calculate the range of cell values in the image for a specific band (start counting from 0).
Definition: ImgReaderGdal.cc:355
ImgReaderGdal::getHistogram
double getHistogram(std::vector< double > &histvector, double &min, double &max, unsigned int &nbin, int theBand=0, bool kde=false)
Calculate the image histogram for a specific band using a defined number of bins and constrained by a...
Definition: ImgReaderGdal.cc:257
ImgReaderGdal::getMin
double getMin(int &col, int &row, int band=0)
Get the minimum cell values for a specific band and report the column and row in which the minimum va...
Definition: ImgReaderGdal.cc:89