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 
28 
30 
36 void ImgReaderGdal::open(const std::string& filename, const GDALAccess& readMode)
37 {
38  m_filename = filename;
39  setCodec(readMode);
40 }
41 
43 {
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 }
~ImgReaderGdal(void)
destructor
int m_nband
number of bands in this dataset
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...
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
bool isNoData(double value) const
Check if value is nodata in this dataset.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
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...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
std::string m_filename
filename of this dataset
std::vector< double > m_noDataValues
no data values for this dataset
double m_gt[6]
geotransform information of this dataset
int nrOfRow(void) const
Get the number of rows of this dataset.
GDALDataset * m_gds
instance of the GDAL dataset of this dataset
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
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 getDeltaY(void) const
Get the pixel cell spacing in y.
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)...
int m_nrow
number of rows in this dataset
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
int m_ncol
number of columns in this dataset
void setCodec(const GDALAccess &readMode=GA_ReadOnly)
Set GDAL dataset number of columns, rows, bands and geotransform.
std::vector< double > m_scale
Vector containing the scale factor to be applied (one scale value for each band)
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...
virtual void close(void)
Close the image.
ImgReaderGdal(void)
default constructor. Image needs to be opened later with one of the open methods. ...
double getDeltaX(void) const
Get the pixel cell spacing in x.
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 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...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98