pktools 2.6.7
Processing Kernel for geospatial data
ImgReaderGdal.cc
1/**********************************************************************
2ImgReaderGdal.cc: class to read raster files using GDAL API library
3Copyright (C) 2008-2016 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 <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
28
30
36void ImgReaderGdal::open(const std::string& filename, const GDALAccess& readMode)
37{
38 m_filename = filename;
39 setCodec(readMode);
40}
41
43{
45}
46
50void 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
89double 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
125double 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
161void 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
208void 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
257double 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
355void 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
373unsigned 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
398void 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}
double getDeltaY(void) const
Get the pixel cell spacing in y.
virtual void close(void)
Close the image.
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) ?
int m_nband
number of bands in this dataset
int m_ncol
number of columns in this dataset
bool isNoData(double value) const
Check if value is nodata in this dataset.
std::string m_filename
filename of this dataset
GDALDataset * m_gds
instance of the GDAL dataset of this dataset
int m_nrow
number of rows in 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::vector< double > m_noDataValues
no data values for this dataset
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
std::vector< double > m_scale
Vector containing the scale factor to be applied (one scale value for each band)
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
double m_gt[6]
geotransform information of this dataset
void setCodec(const GDALAccess &readMode=GA_ReadOnly)
Set GDAL dataset number of columns, rows, bands and geotransform.
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
~ImgReaderGdal(void)
destructor
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 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...
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...
ImgReaderGdal(void)
default constructor. Image needs to be opened later with one of the open methods.
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).
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...
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...