Processing Kernel for remote sensing data
ImgReaderGdal.h
1 /**********************************************************************
2 ImgReaderGdal.h: class to read raster files using GDAL API library
3 Copyright (C) 2008-2012 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 #ifndef _IMGREADERGDAL_H_
21 #define _IMGREADERGDAL_H_
22 
23 #include <assert.h>
24 #include <fstream>
25 #include <string>
26 #include <sstream>
27 #include "gdal_priv.h"
28 #include "base/Vector2d.h"
29 
30 enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };
31 
32 //--------------------------------------------------------------------------
33 class ImgReaderGdal
34 {
35 public:
36  ImgReaderGdal(void);
37  ImgReaderGdal(const std::string& filename){open(filename);};
38  ~ImgReaderGdal(void);
39  void open(const std::string& filename);//, double magicX=1, double magicY=1);
40  void close(void);
41  std::string getFileName() const {return m_filename;};
42  int nrOfCol(void) const { return m_ncol;};
43  int nrOfRow(void) const { return m_nrow;};
44  int nrOfBand(void) const { return m_nband;};
45  bool isGeoRef() const {double gt[6];getGeoTransform(gt);if(gt[5]<0) return true;else return false;};
46  std::string getProjection(void) const;
47  std::string getProjectionRef(void) const;
48  std::string getGeoTransform() const;
49  void getGeoTransform(double* gt) const;
50  /* void getGeoTransform(double& ulx, double& uly, double& deltaX, double& deltaY, double& rot1, double& rot2) const; */
51  std::string getDescription() const;
52  std::string getMetadataItem() const;
53  std::string getImageDescription() const;
54  bool getBoundingBox (double& ulx, double& uly, double& lrx, double& lry) const;
55  bool getCenterPos(double& x, double& y) const;
56  double getUlx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(ulx);};
57  double getUly() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(uly);};
58  double getLrx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lrx);};
59  double getLry() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lry);};
60  // bool getMagicPixel(double& magicX, double& magicY) const {magicX=m_magic_x;magicY=m_magic_y;};
61  void setScale(double theScale, int band=0){
62  /* if(getRasterBand(band)->SetScale(theScale)==CE_Failure){ */
63  if(m_scale.size()!=nrOfBand()){//initialize
64  m_scale.resize(nrOfBand());
65  for(int iband=0;iband<nrOfBand();++iband)
66  m_scale[iband]=1.0;
67  }
68  m_scale[band]=theScale;
69  /* }; */
70  }
71  void setOffset(double theOffset, int band=0){
72  /* if(getRasterBand(band)->SetOffset(theOffset)==CE_Failure){ */
73  if(m_offset.size()!=nrOfBand()){
74  m_offset.resize(nrOfBand());
75  for(int iband=0;iband<nrOfBand();++iband)
76  m_offset[iband]=0.0;
77  }
78  m_offset[band]=theOffset;
79  /* }; */
80  }
81  int getNoDataValues(std::vector<double>& noDataValues) const;
82  bool isNoData(double value) const{if(m_noDataValues.empty()) return false;else return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();};
83  int pushNoDataValue(double noDataValue);
84  int setNoData(const std::vector<double> nodata){m_noDataValues=nodata;};
85  CPLErr GDALSetNoDataValue(double noDataValue, int band=0) {return getRasterBand(band)->SetNoDataValue(noDataValue);};
86  bool covers(double x, double y) const;
87  bool covers(double ulx, double uly, double lrx, double lry) const;
88  bool geo2image(double x, double y, double& i, double& j) const;
89  bool image2geo(double i, double j, double& x, double& y) const;
90  double getDeltaX(void) const {double gt[6];getGeoTransform(gt);return gt[1];};
91  double getDeltaY(void) const {double gt[6];getGeoTransform(gt);return -gt[5];};
92  template<typename T> void readData(T& value, const GDALDataType& dataType, int col, int row, int band=0) const;
93  template<typename T> void readData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int row, int band=0) const;
94  template<typename T> void readData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, double row, int band=0, RESAMPLE resample=NEAR) const;
95  template<typename T> void readDataBlock(Vector2d<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band=0) const;
96  template<typename T> void readDataBlock(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band=0) const;
97  template<typename T> void readData(std::vector<T>& buffer, const GDALDataType& dataType, int row, int band=0) const;
98  template<typename T> void readData(std::vector<T>& buffer, const GDALDataType& dataType, double row, int band=0, RESAMPLE resample=NEAR) const;
99  void getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double& minValue, double& maxValue) const;
100  void getMinMax(double& minValue, double& maxValue, int band=0, bool exhaustiveSearch=false) const;
101  double getMin(int& col, int& row, int band=0) const;
102  unsigned long int getHistogram(std::vector<unsigned long int>& histvector, double& min, double& max,unsigned int& nbin, int theBand=0) const;
103  double getMax(int& col, int& row, int band=0) const;
104  void getRefPix(double& refX, double &refY, int band=0) const;
105  void getRange(std::vector<short>& range, int Band=0) const;
106  GDALDataType getDataType(int band=0) const;
107  GDALRasterBand* getRasterBand(int band=0);
108  GDALColorTable* getColorTable(int band=0) const;
109  std::string getDriverDescription() const;
110  std::string getImageType() const{return getDriverDescription();};
111 // std::string getImageType() const{return "GTiff";};
112  std::string getInterleave() const;
113  std::string getCompression() const;
114  GDALDataset* getDataset(){return m_gds;};
115  char** getMetadata();
116  char** getMetadata() const;
117  void getMetadata(std::list<std::string>& metadata) const;
118 
119 protected:
120  void setCodec();//double magicX, double magicY);
121 
122  std::string m_filename;
123  GDALDataset *m_gds;
124  int m_ncol;
125  int m_nrow;
126  int m_nband;
127  double m_gt[6];
128  /* double m_ulx; */
129  /* double m_uly; */
130  /* double m_delta_x; */
131  /* double m_delta_y; */
132  /* bool m_isGeoRef; */
133  std::vector<double> m_noDataValues;
134  std::vector<double> m_scale;
135  std::vector<double> m_offset;
136 };
137 
138 // adfGeoTransform[0] /* top left x */
139 // adfGeoTransform[1] /* w-e pixel resolution */
140 // adfGeoTransform[2] /* rotation, 0 if image is "north up" */
141 // adfGeoTransform[3] /* top left y */
142 // adfGeoTransform[4] /* rotation, 0 if image is "north up" */
143 // adfGeoTransform[5] /* n-s pixel resolution */
144 
145 template<typename T> void ImgReaderGdal::readData(T& value, const GDALDataType& dataType, int col, int row, int band) const
146 {
147  //fetch raster band
148  GDALRasterBand *poBand;
149  assert(band<nrOfBand()+1);
150  poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
151  assert(col<nrOfCol());
152  assert(col>=0);
153  assert(row<nrOfRow());
154  assert(row>=0);
155  poBand->RasterIO(GF_Read,col,row,1,1,&value,1,1,dataType,0,0);
156  if(m_scale.size()>band)
157  value=static_cast<double>(value)*m_scale[band];
158  if(m_offset.size()>band)
159  value=static_cast<double>(value)+m_offset[band];
160 }
161 
162 template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const GDALDataType& dataType, int minCol, int maxCol, int row, int band) const
163 {
164  //fetch raster band
165  GDALRasterBand *poBand;
166  assert(band<nrOfBand()+1);
167  poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
168  assert(minCol<nrOfCol());
169  assert(minCol>=0);
170  assert(maxCol<nrOfCol());
171  assert(minCol<=maxCol);
172  assert(row<nrOfRow());
173  assert(row>=0);
174  if(buffer.size()!=maxCol-minCol+1)
175  buffer.resize(maxCol-minCol+1);
176  poBand->RasterIO(GF_Read,minCol,row,buffer.size(),1,&(buffer[0]),buffer.size(),1,dataType,0,0);
177  if(m_scale.size()>band||m_offset.size()>band){
178  double theScale=1;
179  double theOffset=0;
180  if(m_scale.size()>band)
181  theScale=m_scale[band];
182  if(m_offset.size()>band)
183  theOffset=m_offset[band];
184  for(int index=0;index<buffer.size();++index)
185  buffer[index]=theScale*static_cast<double>(buffer[index])+theOffset;
186  }
187 }
188 
189 template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, double row, int band, RESAMPLE resample) const
190 {
191  //todo: make upper and lower row depend on isGeo...
192  std::vector<T> readBuffer_upper;
193  std::vector<T> readBuffer_lower;
194  if(buffer.size()!=maxCol-minCol+1)
195  buffer.resize(maxCol-minCol+1);
196  double upperRow=row-0.5;
197  upperRow=static_cast<int>(upperRow);
198  double lowerRow=row+0.5;
199  lowerRow=static_cast<int>(lowerRow);
200  switch(resample){
201  case(BILINEAR):
202  if(lowerRow>=nrOfRow())
203  lowerRow=nrOfRow()-1;
204  if(upperRow<0)
205  upperRow=0;
206  readData(readBuffer_upper,GDT_Float64,minCol,maxCol,static_cast<int>(upperRow),band);
207  readData(readBuffer_lower,GDT_Float64,minCol,maxCol,static_cast<int>(lowerRow),band);
208  //do interpolation in y
209  for(int icol=0;icol<maxCol-minCol+1;++icol){
210  buffer[icol]=(lowerRow-row+0.5)*readBuffer_upper[icol]+(1-lowerRow+row-0.5)*readBuffer_lower[icol];
211  }
212  break;
213  default:
214  readData(buffer,GDT_Float64,minCol,maxCol,static_cast<int>(row),band);
215  break;
216  }
217 }
218 
219 template<typename T> void ImgReaderGdal::readDataBlock(Vector2d<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band) const
220 {
221  buffer.resize(maxRow-minRow+1);
222  for(int irow=minRow;irow<=maxRow;++irow){
223  buffer[irow-minRow].resize(maxCol-minCol+1);
224  readData(buffer[irow-minRow],dataType,minCol,maxCol,irow,band);
225  }
226 }
227 
228 template<typename T> void ImgReaderGdal::readDataBlock(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band) const
229 {
230  double theScale=1;
231  double theOffset=0;
232  if(m_scale.size()>band)
233  theScale=m_scale[band];
234  if(m_offset.size()>band)
235  theOffset=m_offset[band];
236  //fetch raster band
237  GDALRasterBand *poBand;
238  assert(band<nrOfBand()+1);
239  poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
240  if(minCol>=nrOfCol() ||
241  (minCol<0) ||
242  (maxCol>=nrOfCol()) ||
243  (minCol>maxCol) ||
244  (minRow>=nrOfRow()) ||
245  (minRow<0) ||
246  (maxRow>=nrOfRow()) ||
247  (minRow>maxRow)){
248  std::string errorString="block not within image boundaries";
249  throw(errorString);
250  }
251  /* assert(minCol<nrOfCol()); */
252  /* assert(minCol>=0); */
253  /* assert(maxCol<nrOfCol()); */
254  /* assert(minCol<=maxCol); */
255  /* assert(minRow<nrOfRow()); */
256  /* assert(minRow>=0); */
257  /* assert(maxRow<nrOfRow()); */
258  /* assert(minRow<=maxRow); */
259  if(buffer.size()!=(maxRow-minRow+1)*(maxCol-minCol+1))
260  buffer.resize((maxRow-minRow+1)*(maxCol-minCol+1));
261  poBand->RasterIO(GF_Read,minCol,minRow,maxCol-minCol+1,maxRow-minRow+1,&(buffer[0]),(maxCol-minCol+1),(maxRow-minRow+1),dataType,0,0);
262  if(m_scale.size()>band||m_offset.size()>band){
263  for(int index=0;index<buffer.size();++index)
264  buffer[index]=theScale*buffer[index]+theOffset;
265  }
266 }
267 
268 // template<typename T> void ImgReaderGdal::readDataBlock(vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band) const
269 // {
270 // assert(band<nrOfBand()+1);
271 // assert(minCol<nrOfCol());
272 // assert(minCol>=0);
273 // assert(maxCol<nrOfCol());
274 // assert(minCol<=maxCol);
275 // assert(minRow<nrOfRow());
276 // assert(minRow>=0);
277 // assert(maxRow<nrOfRow());
278 // assert(minRow<=maxRow);
279 // if(buffer.size()!=(maxRow-minRow+1)*(maxCol-minCol+1))
280 // buffer.resize((maxRow-minRow+1)*(maxCol-minCol+1));
281 // //fetch raster band
282 // GDALRasterBand *poBand;
283 // assert(band<nrOfBand()+1);
284 // poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
285 // for(int irow=0;irow<maxRow-minRow+1;++irow)
286 // poBand->RasterIO(GF_Read,minCol,minRow+irow,maxCol-minCol+1,1,&(buffer[irow*(maxCol-minCol+1)]),maxCol-minCol+1,1,dataType,0,0);
287 // }
288 
289 template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const GDALDataType& dataType, int row, int band) const
290 {
291  readData(buffer,dataType,0,nrOfCol()-1,row,band);
292 }
293 
294 template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const GDALDataType& dataType, double row, int band, RESAMPLE resample) const
295 {
296  readData(buffer,dataType,0,nrOfCol()-1,row,band,resample);
297 }
298 
299 
300 #endif // _IMGREADERGDAL_H_
301 
302 // //fetch raster band
303 // GDALRasterBand *poBand;
304 // assert(band<nrOfBand()+1);
305 // poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
306 // buffer.resize(maxCol-minCol+1);
307 // assert(minCol<nrOfCol());
308 // assert(row<nrOfRow());
309 // poBand->RasterIO(GF_Read,minCol,row,buffer.size(),1,&(buffer[0]),buffer.size(),1,GDT_Int16,0,0);
Generated on Wed Nov 19 2014 12:26:59 by   doxygen 1.8.6