pktools  2.6.7
Processing Kernel for geospatial data
pkstat.cc
1 /**********************************************************************
2 pkstat.cc: program to calculate basic statistics from raster dataset
3 Copyright (C) 2008-2015 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 <iostream>
21 #include <fstream>
22 #include <math.h>
23 #include "base/Optionpk.h"
24 #include "algorithms/StatFactory.h"
25 #include "algorithms/ImgRegression.h"
26 /******************************************************************************/
78 using namespace std;
79 
80 int main(int argc, char *argv[])
81 {
82  Optionpk<string> input_opt("i","input","name of the input raster dataset");
83  Optionpk<unsigned short> band_opt("b","band","band(s) on which to calculate statistics",0);
84  Optionpk<bool> filename_opt("f", "filename", "Shows image filename ", false);
85  Optionpk<bool> stat_opt("stats", "statistics", "Shows basic statistics (calculate in memory) (min,max, mean and stdDev of the raster datasets)", false);
86  Optionpk<bool> fstat_opt("fstats", "fstatistics", "Shows basic statistics using GDAL computeStatistics (min,max, mean and stdDev of the raster datasets)", false);
87  Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box");
88  Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box");
89  Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box");
90  Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box");
91  Optionpk<double> nodata_opt("nodata","nodata","Set nodata value(s)");
92  Optionpk<short> down_opt("down", "down", "Down sampling factor (for raster sample datasets only). Can be used to create grid points", 1);
93  Optionpk<unsigned int> random_opt("rnd", "rnd", "generate random numbers", 0);
94  Optionpk<double> scale_opt("scale", "scale", "Scale(s) for reading input image(s)");
95  Optionpk<double> offset_opt("offset", "offset", "Offset(s) for reading input image(s)");
96 
97  // Optionpk<bool> transpose_opt("t","transpose","transpose output",false);
98  // Optionpk<std::string> randdist_opt("dist", "dist", "distribution for generating random numbers, see http://www.gn/software/gsl/manual/gsl-ref_toc.html#TOC320 (only uniform and Gaussian supported yet)", "gaussian");
99  // Optionpk<double> randa_opt("rnda", "rnda", "first parameter for random distribution (mean value in case of Gaussian)", 0);
100  // Optionpk<double> randb_opt("rndb", "rndb", "second parameter for random distribution (standard deviation in case of Gaussian)", 1);
101  Optionpk<bool> mean_opt("mean","mean","calculate mean",false);
102  Optionpk<bool> median_opt("median","median","calculate median",false);
103  Optionpk<bool> var_opt("var","var","calculate variance",false);
104  Optionpk<bool> skewness_opt("skew","skewness","calculate skewness",false);
105  Optionpk<bool> kurtosis_opt("kurt","kurtosis","calculate kurtosis",false);
106  Optionpk<bool> stdev_opt("stdev","stdev","calculate standard deviation",false);
107  Optionpk<bool> sum_opt("sum","sum","calculate sum of column",false);
108  Optionpk<bool> minmax_opt("mm","minmax","calculate minimum and maximum value",false);
109  Optionpk<bool> min_opt("min","min","calculate minimum value",false);
110  Optionpk<bool> max_opt("max","max","calculate maximum value",false);
111  Optionpk<double> src_min_opt("src_min","src_min","start reading source from this minimum value");
112  Optionpk<double> src_max_opt("src_max","src_max","stop reading source from this maximum value");
113  Optionpk<bool> histogram_opt("hist","hist","calculate histogram",false);
114  Optionpk<bool> histogram2d_opt("hist2d","hist2d","calculate 2-dimensional histogram based on two images",false);
115  Optionpk<short> nbin_opt("nbin","nbin","number of bins to calculate histogram");
116  Optionpk<bool> relative_opt("rel","relative","use percentiles for histogram to calculate histogram",false);
117  Optionpk<bool> kde_opt("kde","kde","Use Kernel density estimation when producing histogram. The standard deviation is estimated based on Silverman's rule of thumb",false);
118  Optionpk<bool> rmse_opt("rmse","rmse","calculate root mean square error between two raster datasets",false);
119  Optionpk<bool> reg_opt("reg","regression","calculate linear regression between two raster datasets and get correlation coefficient",false);
120  Optionpk<bool> regerr_opt("regerr","regerr","calculate linear regression between two raster datasets and get root mean square error",false);
121  Optionpk<bool> preg_opt("preg","preg","calculate perpendicular regression between two raster datasets and get correlation coefficient",false);
122  Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0,2);
123  fstat_opt.setHide(1);
124  ulx_opt.setHide(1);
125  uly_opt.setHide(1);
126  lrx_opt.setHide(1);
127  lry_opt.setHide(1);
128  down_opt.setHide(1);
129  random_opt.setHide(1);
130  scale_opt.setHide(1);
131  offset_opt.setHide(1);
132  src_min_opt.setHide(1);
133  src_max_opt.setHide(1);
134  kde_opt.setHide(1);
135 
136  // range_opt.setHide(1);
137  // transpose_opt.setHide(1);
138 
139  bool doProcess;//stop process when program was invoked with help option (-h --help)
140  try{
141  //mandatory options
142  doProcess=input_opt.retrieveOption(argc,argv);
143  //optional options
144  band_opt.retrieveOption(argc,argv);
145  filename_opt.retrieveOption(argc,argv);
146  stat_opt.retrieveOption(argc,argv);
147  fstat_opt.retrieveOption(argc,argv);
148  nodata_opt.retrieveOption(argc,argv);
149  mean_opt.retrieveOption(argc,argv);
150  median_opt.retrieveOption(argc,argv);
151  var_opt.retrieveOption(argc,argv);
152  stdev_opt.retrieveOption(argc,argv);
153  minmax_opt.retrieveOption(argc,argv);
154  min_opt.retrieveOption(argc,argv);
155  max_opt.retrieveOption(argc,argv);
156  histogram_opt.retrieveOption(argc,argv);
157  nbin_opt.retrieveOption(argc,argv);
158  relative_opt.retrieveOption(argc,argv);
159  histogram2d_opt.retrieveOption(argc,argv);
160  rmse_opt.retrieveOption(argc,argv);
161  reg_opt.retrieveOption(argc,argv);
162  regerr_opt.retrieveOption(argc,argv);
163  preg_opt.retrieveOption(argc,argv);
164  //advanced options
165  ulx_opt.retrieveOption(argc,argv);
166  uly_opt.retrieveOption(argc,argv);
167  lrx_opt.retrieveOption(argc,argv);
168  lry_opt.retrieveOption(argc,argv);
169  down_opt.retrieveOption(argc,argv);
170  random_opt.retrieveOption(argc,argv);
171  scale_opt.retrieveOption(argc,argv);
172  offset_opt.retrieveOption(argc,argv);
173  src_min_opt.retrieveOption(argc,argv);
174  src_max_opt.retrieveOption(argc,argv);
175  kde_opt.retrieveOption(argc,argv);
176  verbose_opt.retrieveOption(argc,argv);
177  }
178  catch(string predefinedString){
179  std::cout << predefinedString << std::endl;
180  exit(0);
181  }
182  if(!doProcess){
183  cout << endl;
184  cout << "Usage: pkstat -i input" << endl;
185  cout << endl;
186  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
187  exit(0);//help was invoked, stop processing
188  }
189 
190  if(src_min_opt.size()){
191  while(src_min_opt.size()<band_opt.size())
192  src_min_opt.push_back(src_min_opt[0]);
193  }
194  if(src_max_opt.size()){
195  while(src_max_opt.size()<band_opt.size())
196  src_max_opt.push_back(src_max_opt[0]);
197  }
198 
199  unsigned int nbin=0;
200  double minX=0;
201  double minY=0;
202  double maxX=0;
203  double maxY=0;
204  double minValue=(src_min_opt.size())? src_min_opt[0] : 0;
205  double maxValue=(src_max_opt.size())? src_max_opt[0] : 0;
206  double meanValue=0;
207  double medianValue=0;
208  double stdDev=0;
209 
210  const char* pszMessage;
211  void* pProgressArg=NULL;
212  GDALProgressFunc pfnProgress=GDALTermProgress;
213  double progress=0;
214  srand(time(NULL));
215 
218  std::vector<double> histogramOutput;
219  double nsample=0;
220 
221  ImgReaderGdal imgReader;
222 
223  if(scale_opt.size()){
224  while(scale_opt.size()<input_opt.size())
225  scale_opt.push_back(scale_opt[0]);
226  }
227  if(offset_opt.size()){
228  while(offset_opt.size()<input_opt.size())
229  offset_opt.push_back(offset_opt[0]);
230  }
231  if(input_opt.empty()){
232  std::cerr << "No image dataset provided (use option -i). Use --help for help information";
233  exit(0);
234  }
235  for(int ifile=0;ifile<input_opt.size();++ifile){
236  try{
237  imgReader.open(input_opt[ifile]);
238  }
239  catch(std::string errorstring){
240  std::cout << errorstring << std::endl;
241  exit(0);
242  }
243 
244  if(filename_opt[0])
245  std::cout << " --input " << input_opt[ifile] << " ";
246 
247  for(int inodata=0;inodata<nodata_opt.size();++inodata)
248  imgReader.pushNoDataValue(nodata_opt[inodata]);
249 
250  int nband=band_opt.size();
251  for(int iband=0;iband<nband;++iband){
252  minValue=(src_min_opt.size())? src_min_opt[0] : 0;
253  maxValue=(src_max_opt.size())? src_max_opt[0] : 0;
254  for(int inodata=0;inodata<nodata_opt.size();++inodata){
255  if(!inodata)
256  imgReader.GDALSetNoDataValue(nodata_opt[0],band_opt[iband]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
257  }
258 
259  if(offset_opt.size()>ifile)
260  imgReader.setOffset(offset_opt[ifile],band_opt[iband]);
261  if(scale_opt.size()>ifile)
262  imgReader.setScale(scale_opt[ifile],band_opt[iband]);
263 
264  if(stat_opt[0]||mean_opt[0]||median_opt[0]||var_opt[0]||stdev_opt[0]){//the hard way (in memory)
266  vector<double> readBuffer;
267  double varValue;
268  imgReader.readDataBlock(readBuffer, 0, imgReader.nrOfCol()-1, 0, imgReader.nrOfRow()-1, band_opt[0]);
269  stat.setNoDataValues(nodata_opt);
270  stat.meanVar(readBuffer,meanValue,varValue);
271  medianValue=stat.median(readBuffer);
272  stat.minmax(readBuffer,readBuffer.begin(),readBuffer.end(),minValue,maxValue);
273  if(mean_opt[0])
274  std::cout << "--mean " << meanValue << " ";
275  if(median_opt[0])
276  std::cout << "--median " << medianValue << " ";
277  if(stdev_opt[0])
278  std::cout << "--stdDev " << sqrt(varValue) << " ";
279  if(var_opt[0])
280  std::cout << "--var " << varValue << " ";
281  if(stat_opt[0])
282  std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << sqrt(varValue) << " ";
283  }
284 
285  if(fstat_opt[0]){//the fast way
286  assert(band_opt[iband]<imgReader.nrOfBand());
287  GDALProgressFunc pfnProgress;
288  void* pProgressData;
289  GDALRasterBand* rasterBand;
290  rasterBand=imgReader.getRasterBand(band_opt[iband]);
291  rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
292 
293  std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " ";
294  }
295 
296  if(minmax_opt[0]||min_opt[0]||max_opt[0]){
297  assert(band_opt[iband]<imgReader.nrOfBand());
298 
299  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]))){
300  double uli,ulj,lri,lrj;
301  imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
302  imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
303  imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue);
304  }
305  else{
306  imgReader.getMinMax(minValue,maxValue,band_opt[iband]);
307  }
308  if(minmax_opt[0])
309  std::cout << "-min " << minValue << " -max " << maxValue << " ";
310  else{
311  if(min_opt[0])
312  std::cout << "-min " << minValue << " ";
313  if(max_opt[0])
314  std::cout << "-max " << maxValue << " ";
315  }
316  }
317  }
318  if(histogram_opt[0]){//aggregate results from multiple inputs, but only calculate for first selected band
319  assert(band_opt[0]<imgReader.nrOfBand());
320  nbin=(nbin_opt.size())? nbin_opt[0]:0;
321 
322  imgReader.getMinMax(minValue,maxValue,band_opt[0]);
323  if(src_min_opt.size())
324  minValue=src_min_opt[0];
325  if(src_max_opt.size())
326  maxValue=src_max_opt[0];
327  if(minValue>=maxValue)
328  imgReader.getMinMax(minValue,maxValue,band_opt[0]);
329 
330  if(verbose_opt[0])
331  cout << "number of valid pixels in image: " << imgReader.getNvalid(band_opt[0]) << endl;
332 
333  nsample+=imgReader.getHistogram(histogramOutput,minValue,maxValue,nbin,band_opt[0],kde_opt[0]);
334 
335  //only output for last input file
336  if(ifile==input_opt.size()-1){
337  std::cout.precision(10);
338  for(int bin=0;bin<nbin;++bin){
339  double binValue=0;
340  if(nbin==maxValue-minValue+1)
341  binValue=minValue+bin;
342  else
343  binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin;
344  std::cout << binValue << " ";
345  if(relative_opt[0]||kde_opt[0])
346  std::cout << 100.0*static_cast<double>(histogramOutput[bin])/static_cast<double>(nsample) << std::endl;
347  else
348  std::cout << static_cast<double>(histogramOutput[bin]) << std::endl;
349  }
350  }
351  }
352  if(histogram2d_opt[0]&&input_opt.size()<2){
353  assert(band_opt.size()>1);
354  imgReader.getMinMax(minX,maxX,band_opt[0]);
355  imgReader.getMinMax(minY,maxY,band_opt[1]);
356  if(src_min_opt.size()){
357  minX=src_min_opt[0];
358  minY=src_min_opt[1];
359  }
360  if(src_max_opt.size()){
361  maxX=src_max_opt[0];
362  maxY=src_max_opt[1];
363  }
364  nbin=(nbin_opt.size())? nbin_opt[0]:0;
365  if(nbin<=1){
366  std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
367  if(minX>=maxX)
368  imgReader.getMinMax(minX,maxX,band_opt[0]);
369  if(minY>=maxY)
370  imgReader.getMinMax(minY,maxY,band_opt[1]);
371 
372  minValue=(minX<minY)? minX:minY;
373  maxValue=(maxX>maxY)? maxX:maxY;
374  if(verbose_opt[0])
375  std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
376  nbin=maxValue-minValue+1;
377  }
378  assert(nbin>1);
379  double sigma=0;
380  //kernel density estimation as in http://en.wikipedia.org/wiki/Kernel_density_estimation
381  if(kde_opt[0]){
382  assert(band_opt[0]<imgReader.nrOfBand());
383  assert(band_opt[1]<imgReader.nrOfBand());
384  GDALProgressFunc pfnProgress;
385  void* pProgressData;
386  GDALRasterBand* rasterBand;
387  double stdDev1=0;
388  double stdDev2=0;
389  rasterBand=imgReader.getRasterBand(band_opt[0]);
390  rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData);
391  rasterBand=imgReader.getRasterBand(band_opt[1]);
392  rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData);
393 
394  double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
395  if(random_opt[0]>0)
396  estimatedSize*=random_opt[0]/100.0;
397  sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2);
398  }
399  assert(nbin);
400  if(verbose_opt[0]){
401  if(sigma>0)
402  std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for bands " << band_opt[0] << " and " << band_opt[1] << std::endl;
403  else
404  std::cout << "calculating 2d histogram for bands " << band_opt[0] << " and " << band_opt[1] << std::endl;
405  std::cout << "nbin: " << nbin << std::endl;
406  }
407 
408 
409  vector< vector<double> > output;
410 
411  if(maxX<=minX)
412  imgReader.getMinMax(minX,maxX,band_opt[0]);
413  if(maxY<=minY)
414  imgReader.getMinMax(minY,maxY,band_opt[1]);
415 
416  if(maxX<=minX){
417  std::ostringstream s;
418  s<<"Error: could not calculate distribution (minX>=maxX)";
419  throw(s.str());
420  }
421  if(maxY<=minY){
422  std::ostringstream s;
423  s<<"Error: could not calculate distribution (minY>=maxY)";
424  throw(s.str());
425  }
426  output.resize(nbin);
427  for(int i=0;i<nbin;++i){
428  output[i].resize(nbin);
429  for(int j=0;j<nbin;++j)
430  output[i][j]=0;
431  }
432  int binX=0;
433  int binY=0;
434  vector<double> inputX(imgReader.nrOfCol());
435  vector<double> inputY(imgReader.nrOfCol());
436  unsigned long int nvalid=0;
437  for(int irow=0;irow<imgReader.nrOfRow();++irow){
438  if(irow%down_opt[0])
439  continue;
440  imgReader.readData(inputX,irow,band_opt[0]);
441  imgReader.readData(inputY,irow,band_opt[1]);
442  for(int icol=0;icol<imgReader.nrOfCol();++icol){
443  if(icol%down_opt[0])
444  continue;
445  if(random_opt[0]>0){
446  double p=static_cast<double>(rand())/(RAND_MAX);
447  p*=100.0;
448  if(p>random_opt[0])
449  continue;//do not select for now, go to next column
450  }
451  if(imgReader.isNoData(inputX[icol]))
452  continue;
453  if(imgReader.isNoData(inputY[icol]))
454  continue;
455  ++nvalid;
456  if(inputX[icol]>=maxX)
457  binX=nbin-1;
458  else if(inputX[icol]<=minX)
459  binX=0;
460  else
461  binX=static_cast<int>(static_cast<double>(inputX[icol]-minX)/(maxX-minX)*nbin);
462  if(inputY[icol]>=maxY)
463  binY=nbin-1;
464  else if(inputY[icol]<=minX)
465  binY=0;
466  else
467  binY=static_cast<int>(static_cast<double>(inputY[icol]-minY)/(maxY-minY)*nbin);
468  assert(binX>=0);
469  assert(binX<output.size());
470  assert(binY>=0);
471  assert(binY<output[binX].size());
472  if(sigma>0){
473  //create kde for Gaussian basis function
474  //todo: speed up by calculating first and last bin with non-zero contriubtion...
475  for(int ibinX=0;ibinX<nbin;++ibinX){
476  double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
477  double pdfX=gsl_ran_gaussian_pdf(inputX[icol]-centerX, sigma);
478  for(int ibinY=0;ibinY<nbin;++ibinY){
479  //calculate \integral_ibinX^(ibinX+1)
480  double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
481  double pdfY=gsl_ran_gaussian_pdf(inputY[icol]-centerY, sigma);
482  output[ibinX][binY]+=pdfX*pdfY;
483  }
484  }
485  }
486  else
487  ++output[binX][binY];
488  }
489  }
490  if(verbose_opt[0])
491  cout << "number of valid pixels: " << nvalid << endl;
492 
493  for(int binX=0;binX<nbin;++binX){
494  cout << endl;
495  for(int binY=0;binY<nbin;++binY){
496  double binValueX=0;
497  if(nbin==maxX-minX+1)
498  binValueX=minX+binX;
499  else
500  binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
501  double binValueY=0;
502  if(nbin==maxY-minY+1)
503  binValueY=minY+binY;
504  else
505  binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
506 
507  double value=static_cast<double>(output[binX][binY]);
508 
509  if(relative_opt[0])
510  value*=100.0/nvalid;
511 
512  cout << binValueX << " " << binValueY << " " << value << std::endl;
513  // double value=static_cast<double>(output[binX][binY])/nvalid;
514  // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
515  }
516  }
517  }
518  if(reg_opt[0]&&input_opt.size()<2){
519  if(band_opt.size()<2)
520  continue;
521  imgreg.setDown(down_opt[0]);
522  imgreg.setThreshold(random_opt[0]);
523  double c0=0;//offset
524  double c1=1;//scale
525  double r2=imgreg.getR2(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]);
526  std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl;
527  }
528  if(regerr_opt[0]&&input_opt.size()<2){
529  if(band_opt.size()<2)
530  continue;
531  imgreg.setDown(down_opt[0]);
532  imgreg.setThreshold(random_opt[0]);
533  double c0=0;//offset
534  double c1=1;//scale
535  double err=imgreg.getRMSE(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]);
536  std::cout << "-c0 " << c0 << " -c1 " << c1 << " -rmse " << err << std::endl;
537  }
538  if(rmse_opt[0]&&input_opt.size()<2){
539  if(band_opt.size()<2)
540  continue;
541  vector<double> xBuffer(imgReader.nrOfCol());
542  vector<double> yBuffer(imgReader.nrOfCol());
543  double mse=0;
544  double nValid=0;
545  double nPixel=imgReader.nrOfCol()/down_opt[0]*imgReader.nrOfRow()/down_opt[0];
546  for(int irow;irow<imgReader.nrOfRow();irow+=down_opt[0]){
547  imgReader.readData(xBuffer,irow,band_opt[0]);
548  imgReader.readData(yBuffer,irow,band_opt[1]);
549  for(int icol;icol<imgReader.nrOfCol();icol+=down_opt[0]){
550  double xValue=xBuffer[icol];
551  double yValue=yBuffer[icol];
552  if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
553  continue;
554  }
555  if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
556  continue;
557  }
558  if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[0]||yValue>src_max_opt[0])
559  continue;
560  ++nValid;
561  double e=xValue-yValue;
562  if(relative_opt[0])
563  e/=yValue;
564  mse+=e*e/nPixel;
565  }
566  }
567  double correctNorm=nValid;
568  correctNorm/=nPixel;
569  mse/=correctNorm;
570  std::cout << " -rmse " << sqrt(mse) << std::endl;
571  }
572  if(preg_opt[0]&&input_opt.size()<2){
573  if(band_opt.size()<2)
574  continue;
575  imgreg.setDown(down_opt[0]);
576  imgreg.setThreshold(random_opt[0]);
577  double c0=0;//offset
578  double c1=1;//scale
579  double r2=imgreg.pgetR2(imgReader,band_opt[0],band_opt[1],c0,c1,verbose_opt[0]);
580  std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl;
581  }
582  imgReader.close();
583  }
584  // if(rmse_opt[0]&&(input_opt.size()>1)){
585  // while(band_opt.size()<input_opt.size())
586  // band_opt.push_back(band_opt[0]);
587  // if(src_min_opt.size()){
588  // while(src_min_opt.size()<input_opt.size())
589  // src_min_opt.push_back(src_min_opt[0]);
590  // }
591  // if(src_max_opt.size()){
592  // while(src_max_opt.size()<input_opt.size())
593  // src_max_opt.push_back(src_max_opt[0]);
594  // }
595  // ImgReaderGdal imgReader1(input_opt[0]);
596  // ImgReaderGdal imgReader2(input_opt[1]);
597 
598  // if(offset_opt.size())
599  // imgReader1.setOffset(offset_opt[0],band_opt[0]);
600  // if(scale_opt.size())
601  // imgReader1.setScale(scale_opt[0],band_opt[0]);
602  // if(offset_opt.size()>1)
603  // imgReader2.setOffset(offset_opt[1],band_opt[1]);
604  // if(scale_opt.size()>1)
605  // imgReader2.setScale(scale_opt[1],band_opt[1]);
606 
607  // for(int inodata=0;inodata<nodata_opt.size();++inodata){
608  // imgReader1.pushNoDataValue(nodata_opt[inodata]);
609  // imgReader2.pushNoDataValue(nodata_opt[inodata]);
610  // }
611  // vector<double> xBuffer(imgReader1.nrOfCol());
612  // vector<double> yBuffer(imgReader2.nrOfCol());
613  // double mse=0;
614  // double nValid=0;
615  // double nPixel=imgReader.nrOfCol()/imgReader.nrOfRow()/down_opt[0]/down_opt[0];
616  // for(int irow;irow<imgReader1.nrOfRow();irow+=down_opt[0]){
617  // double irow1=irow;
618  // double irow2=0;
619  // double icol1=0;
620  // double icol2=0;
621  // double geoX=0;
622  // double geoY=0;
623  // imgReader1.image2geo(icol1,irow1,geoX,geoY);
624  // imgReader2.geo2image(geoX,geoY,icol2,irow2);
625  // irow2=static_cast<int>(irow2);
626  // imgReader1.readData(xBuffer,irow1,band_opt[0]);
627  // imgReader2.readData(yBuffer,irow2,band_opt[1]);
628  // for(int icol;icol<imgReader.nrOfCol();icol+=down_opt[0]){
629  // icol1=icol;
630  // imgReader1.image2geo(icol1,irow1,geoX,geoY);
631  // imgReader2.geo2image(geoX,geoY,icol2,irow2);
632  // double xValue=xBuffer[icol1];
633  // double yValue=yBuffer[icol2];
634  // if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
635  // continue;
636  // }
637  // if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[1]||yValue>src_max_opt[1])
638  // continue;
639  // ++nValid;
640  // double e=xValue-yValue;
641  // if(relative_opt[0])
642  // e/=yValue;
643  // mse+=e*e/nPixel;
644  // }
645  // }
646  // double correctNorm=nValid;
647  // correctNorm/=nPixel;
648  // mse/=correctNorm;
649  // std::cout << " -rmse " << sqrt(mse) << std::endl;
650  // }
651  if(reg_opt[0]&&(input_opt.size()>1)){
652  imgreg.setDown(down_opt[0]);
653  imgreg.setThreshold(random_opt[0]);
654  double c0=0;//offset
655  double c1=1;//scale
656  while(band_opt.size()<input_opt.size())
657  band_opt.push_back(band_opt[0]);
658  if(src_min_opt.size()){
659  while(src_min_opt.size()<input_opt.size())
660  src_min_opt.push_back(src_min_opt[0]);
661  }
662  if(src_max_opt.size()){
663  while(src_max_opt.size()<input_opt.size())
664  src_max_opt.push_back(src_max_opt[0]);
665  }
666  ImgReaderGdal imgReader1(input_opt[0]);
667  ImgReaderGdal imgReader2(input_opt[1]);
668 
669  if(offset_opt.size())
670  imgReader1.setOffset(offset_opt[0],band_opt[0]);
671  if(scale_opt.size())
672  imgReader1.setScale(scale_opt[0],band_opt[0]);
673  if(offset_opt.size()>1)
674  imgReader2.setOffset(offset_opt[1],band_opt[1]);
675  if(scale_opt.size()>1)
676  imgReader2.setScale(scale_opt[1],band_opt[1]);
677 
678  for(int inodata=0;inodata<nodata_opt.size();++inodata){
679  if(!inodata){
680  imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
681  imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
682  }
683  imgReader1.pushNoDataValue(nodata_opt[inodata]);
684  imgReader2.pushNoDataValue(nodata_opt[inodata]);
685  }
686 
687  double r2=imgreg.getR2(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]);
688  std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl;
689  imgReader1.close();
690  imgReader2.close();
691  }
692  if(preg_opt[0]&&(input_opt.size()>1)){
693  imgreg.setDown(down_opt[0]);
694  imgreg.setThreshold(random_opt[0]);
695  double c0=0;//offset
696  double c1=1;//scale
697  while(band_opt.size()<input_opt.size())
698  band_opt.push_back(band_opt[0]);
699  if(src_min_opt.size()){
700  while(src_min_opt.size()<input_opt.size())
701  src_min_opt.push_back(src_min_opt[0]);
702  }
703  if(src_max_opt.size()){
704  while(src_max_opt.size()<input_opt.size())
705  src_max_opt.push_back(src_max_opt[0]);
706  }
707  ImgReaderGdal imgReader1(input_opt[0]);
708  ImgReaderGdal imgReader2(input_opt[1]);
709 
710  if(offset_opt.size())
711  imgReader1.setOffset(offset_opt[0],band_opt[0]);
712  if(scale_opt.size())
713  imgReader1.setScale(scale_opt[0],band_opt[0]);
714  if(offset_opt.size()>1)
715  imgReader2.setOffset(offset_opt[1],band_opt[1]);
716  if(scale_opt.size()>1)
717  imgReader2.setScale(scale_opt[1],band_opt[1]);
718 
719  for(int inodata=0;inodata<nodata_opt.size();++inodata){
720  if(!inodata){
721  imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
722  imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
723  }
724  imgReader1.pushNoDataValue(nodata_opt[inodata]);
725  imgReader2.pushNoDataValue(nodata_opt[inodata]);
726  }
727 
728  double r2=imgreg.pgetR2(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]);
729  std::cout << "-c0 " << c0 << " -c1 " << c1 << " -r2 " << r2 << std::endl;
730  imgReader1.close();
731  imgReader2.close();
732  }
733  if(regerr_opt[0]&&(input_opt.size()>1)){
734  imgreg.setDown(down_opt[0]);
735  imgreg.setThreshold(random_opt[0]);
736  double c0=0;//offset
737  double c1=1;//scale
738  while(band_opt.size()<input_opt.size())
739  band_opt.push_back(band_opt[0]);
740  if(src_min_opt.size()){
741  while(src_min_opt.size()<input_opt.size())
742  src_min_opt.push_back(src_min_opt[0]);
743  }
744  if(src_max_opt.size()){
745  while(src_max_opt.size()<input_opt.size())
746  src_max_opt.push_back(src_max_opt[0]);
747  }
748  ImgReaderGdal imgReader1(input_opt[0]);
749  ImgReaderGdal imgReader2(input_opt[1]);
750 
751  if(offset_opt.size())
752  imgReader1.setOffset(offset_opt[0],band_opt[0]);
753  if(scale_opt.size())
754  imgReader1.setScale(scale_opt[0],band_opt[0]);
755  if(offset_opt.size()>1)
756  imgReader2.setOffset(offset_opt[1],band_opt[1]);
757  if(scale_opt.size()>1)
758  imgReader2.setScale(scale_opt[1],band_opt[1]);
759 
760  for(int inodata=0;inodata<nodata_opt.size();++inodata){
761  if(!inodata){
762  imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
763  imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
764  }
765  imgReader1.pushNoDataValue(nodata_opt[inodata]);
766  imgReader2.pushNoDataValue(nodata_opt[inodata]);
767  }
768 
769  double err=imgreg.getRMSE(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]);
770  std::cout << "-c0 " << c0 << " -c1 " << c1 << " -rmse " << err << std::endl;
771  imgReader1.close();
772  imgReader2.close();
773  }
774  if(rmse_opt[0]&&(input_opt.size()>1)){
775  imgreg.setDown(down_opt[0]);
776  imgreg.setThreshold(random_opt[0]);
777  double c0=0;//offset
778  double c1=1;//scale
779  while(band_opt.size()<input_opt.size())
780  band_opt.push_back(band_opt[0]);
781  if(src_min_opt.size()){
782  while(src_min_opt.size()<input_opt.size())
783  src_min_opt.push_back(src_min_opt[0]);
784  }
785  if(src_max_opt.size()){
786  while(src_max_opt.size()<input_opt.size())
787  src_max_opt.push_back(src_max_opt[0]);
788  }
789  ImgReaderGdal imgReader1(input_opt[0]);
790  ImgReaderGdal imgReader2(input_opt[1]);
791 
792  if(offset_opt.size())
793  imgReader1.setOffset(offset_opt[0],band_opt[0]);
794  if(scale_opt.size())
795  imgReader1.setScale(scale_opt[0],band_opt[0]);
796  if(offset_opt.size()>1)
797  imgReader2.setOffset(offset_opt[1],band_opt[1]);
798  if(scale_opt.size()>1)
799  imgReader2.setScale(scale_opt[1],band_opt[1]);
800 
801  for(int inodata=0;inodata<nodata_opt.size();++inodata){
802  if(!inodata){
803  imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
804  imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
805  }
806  imgReader1.pushNoDataValue(nodata_opt[inodata]);
807  imgReader2.pushNoDataValue(nodata_opt[inodata]);
808  }
809 
810  double err=imgreg.getRMSE(imgReader1,imgReader2,c0,c1,band_opt[0],band_opt[1],verbose_opt[0]);
811  std::cout << "-rmse " << err << std::endl;
812  imgReader1.close();
813  imgReader2.close();
814  }
815  if(histogram2d_opt[0]&&(input_opt.size()>1)){
816  while(band_opt.size()<input_opt.size())
817  band_opt.push_back(band_opt[0]);
818  if(src_min_opt.size()){
819  while(src_min_opt.size()<input_opt.size())
820  src_min_opt.push_back(src_min_opt[0]);
821  }
822  if(src_max_opt.size()){
823  while(src_max_opt.size()<input_opt.size())
824  src_max_opt.push_back(src_max_opt[0]);
825  }
826  ImgReaderGdal imgReader1(input_opt[0]);
827  ImgReaderGdal imgReader2(input_opt[1]);
828 
829  if(offset_opt.size())
830  imgReader1.setOffset(offset_opt[0],band_opt[0]);
831  if(scale_opt.size())
832  imgReader1.setScale(scale_opt[0],band_opt[0]);
833  if(offset_opt.size()>1)
834  imgReader2.setOffset(offset_opt[1],band_opt[1]);
835  if(scale_opt.size()>1)
836  imgReader2.setScale(scale_opt[1],band_opt[1]);
837 
838  for(int inodata=0;inodata<nodata_opt.size();++inodata){
839  if(!inodata){
840  imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
841  imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
842  }
843  imgReader1.pushNoDataValue(nodata_opt[inodata]);
844  imgReader2.pushNoDataValue(nodata_opt[inodata]);
845  }
846 
847  imgReader1.getMinMax(minX,maxX,band_opt[0]);
848  imgReader2.getMinMax(minY,maxY,band_opt[1]);
849 
850  if(verbose_opt[0]){
851  cout << "minX: " << minX << endl;
852  cout << "maxX: " << maxX << endl;
853  cout << "minY: " << minY << endl;
854  cout << "maxY: " << maxY << endl;
855  }
856 
857  if(src_min_opt.size()){
858  minX=src_min_opt[0];
859  minY=src_min_opt[1];
860  }
861  if(src_max_opt.size()){
862  maxX=src_max_opt[0];
863  maxY=src_max_opt[1];
864  }
865 
866  nbin=(nbin_opt.size())? nbin_opt[0]:0;
867  if(nbin<=1){
868  std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
869  // imgReader1.getMinMax(minX,maxX,band_opt[0]);
870  // imgReader2.getMinMax(minY,maxY,band_opt[0]);
871  if(minX>=maxX)
872  imgReader1.getMinMax(minX,maxX,band_opt[0]);
873  if(minY>=maxY)
874  imgReader2.getMinMax(minY,maxY,band_opt[1]);
875 
876  minValue=(minX<minY)? minX:minY;
877  maxValue=(maxX>maxY)? maxX:maxY;
878  if(verbose_opt[0])
879  std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
880  nbin=maxValue-minValue+1;
881  }
882  assert(nbin>1);
883  double sigma=0;
884  //kernel density estimation as in http://en.wikipedia.org/wiki/Kernel_density_estimation
885  if(kde_opt[0]){
886  GDALProgressFunc pfnProgress;
887  void* pProgressData;
888  GDALRasterBand* rasterBand;
889  double stdDev1=0;
890  double stdDev2=0;
891  rasterBand=imgReader1.getRasterBand(band_opt[0]);
892  rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData);
893  rasterBand=imgReader2.getRasterBand(band_opt[0]);
894  rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData);
895 
896  //todo: think of smarter way how to estimate size (nodata!)
897  double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
898  if(random_opt[0]>0)
899  estimatedSize*=random_opt[0]/100.0;
900  sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2);
901  }
902  assert(nbin);
903  if(verbose_opt[0]){
904  if(sigma>0)
905  std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl;
906  else
907  std::cout << "calculating 2d histogram for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl;
908  std::cout << "nbin: " << nbin << std::endl;
909  }
910 
911  vector< vector<double> > output;
912 
913  if(maxX<=minX)
914  imgReader1.getMinMax(minX,maxX,band_opt[0]);
915  if(maxY<=minY)
916  imgReader2.getMinMax(minY,maxY,band_opt[1]);
917 
918  if(maxX<=minX){
919  std::ostringstream s;
920  s<<"Error: could not calculate distribution (minX>=maxX)";
921  throw(s.str());
922  }
923  if(maxY<=minY){
924  std::ostringstream s;
925  s<<"Error: could not calculate distribution (minY>=maxY)";
926  throw(s.str());
927  }
928  if(verbose_opt[0]){
929  cout << "minX: " << minX << endl;
930  cout << "maxX: " << maxX << endl;
931  cout << "minY: " << minY << endl;
932  cout << "maxY: " << maxY << endl;
933  }
934  output.resize(nbin);
935  for(int i=0;i<nbin;++i){
936  output[i].resize(nbin);
937  for(int j=0;j<nbin;++j)
938  output[i][j]=0;
939  }
940  int binX=0;
941  int binY=0;
942  vector<double> inputX(imgReader1.nrOfCol());
943  vector<double> inputY(imgReader2.nrOfCol());
944  double nvalid=0;
945  double geoX=0;
946  double geoY=0;
947  double icol1=0;
948  double irow1=0;
949  double icol2=0;
950  double irow2=0;
951  for(int irow=0;irow<imgReader1.nrOfRow();++irow){
952  if(irow%down_opt[0])
953  continue;
954  irow1=irow;
955  imgReader1.image2geo(icol1,irow1,geoX,geoY);
956  imgReader2.geo2image(geoX,geoY,icol2,irow2);
957  irow2=static_cast<int>(irow2);
958  imgReader1.readData(inputX,irow1,band_opt[0]);
959  imgReader2.readData(inputY,irow2,band_opt[1]);
960  for(int icol=0;icol<imgReader.nrOfCol();++icol){
961  if(icol%down_opt[0])
962  continue;
963  icol1=icol;
964  if(random_opt[0]>0){
965  double p=static_cast<double>(rand())/(RAND_MAX);
966  p*=100.0;
967  if(p>random_opt[0])
968  continue;//do not select for now, go to next column
969  }
970  if(imgReader1.isNoData(inputX[icol]))
971  continue;
972  imgReader1.image2geo(icol1,irow1,geoX,geoY);
973  imgReader2.geo2image(geoX,geoY,icol2,irow2);
974  icol2=static_cast<int>(icol2);
975  if(imgReader2.isNoData(inputY[icol2]))
976  continue;
977  // ++nvalid;
978  if(inputX[icol1]>=maxX)
979  binX=nbin-1;
980  else if(inputX[icol]<=minX)
981  binX=0;
982  else
983  binX=static_cast<int>(static_cast<double>(inputX[icol1]-minX)/(maxX-minX)*nbin);
984  if(inputY[icol2]>=maxY)
985  binY=nbin-1;
986  else if(inputY[icol2]<=minY)
987  binY=0;
988  else
989  binY=static_cast<int>(static_cast<double>(inputY[icol2]-minY)/(maxY-minY)*nbin);
990  assert(binX>=0);
991  assert(binX<output.size());
992  assert(binY>=0);
993  assert(binY<output[binX].size());
994  if(sigma>0){
995  //create kde for Gaussian basis function
996  //todo: speed up by calculating first and last bin with non-zero contriubtion...
997  for(int ibinX=0;ibinX<nbin;++ibinX){
998  double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
999  double pdfX=gsl_ran_gaussian_pdf(inputX[icol1]-centerX, sigma);
1000  for(int ibinY=0;ibinY<nbin;++ibinY){
1001  //calculate \integral_ibinX^(ibinX+1)
1002  double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
1003  double pdfY=gsl_ran_gaussian_pdf(inputY[icol2]-centerY, sigma);
1004  output[ibinX][binY]+=pdfX*pdfY;
1005  nvalid+=pdfX*pdfY;
1006  }
1007  }
1008  }
1009  else{
1010  ++output[binX][binY];
1011  ++nvalid;
1012  }
1013  }
1014  }
1015  if(verbose_opt[0])
1016  cout << "number of valid pixels: " << nvalid << endl;
1017  for(int binX=0;binX<nbin;++binX){
1018  cout << endl;
1019  for(int binY=0;binY<nbin;++binY){
1020  double binValueX=0;
1021  if(nbin==maxX-minX+1)
1022  binValueX=minX+binX;
1023  else
1024  binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
1025  double binValueY=0;
1026  if(nbin==maxY-minY+1)
1027  binValueY=minY+binY;
1028  else
1029  binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
1030  double value=static_cast<double>(output[binX][binY]);
1031 
1032  if(relative_opt[0]||kde_opt[0])
1033  value*=100.0/nvalid;
1034 
1035  cout << binValueX << " " << binValueY << " " << value << std::endl;
1036  // double value=static_cast<double>(output[binX][binY])/nvalid;
1037  // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
1038  }
1039  }
1040  imgReader1.close();
1041  imgReader2.close();
1042  }
1043 
1044  if(!histogram_opt[0]||histogram2d_opt[0])
1045  std::cout << std::endl;
1046 }
1047 
1048 // int nband=(band_opt.size()) ? band_opt.size() : imgReader.nrOfBand();
1049 
1050 // const char* pszMessage;
1051 // void* pProgressArg=NULL;
1052 // GDALProgressFunc pfnProgress=GDALTermProgress;
1053 // double progress=0;
1054 // srand(time(NULL));
1055 
1056 
1057 // statfactory::StatFactory stat;
1058 // imgregression::ImgRegression imgreg;
1059 
1060 // pfnProgress(progress,pszMessage,pProgressArg);
1061 // for(irow=0;irow<classReader.nrOfRow();++irow){
1062 // if(irow%down_opt[0])
1063 // continue;
1064 // // classReader.readData(classBuffer,irow);
1065 // classReader.readData(classBuffer,irow);
1066 // double x,y;//geo coordinates
1067 // double iimg,jimg;//image coordinates in img image
1068 // for(icol=0;icol<classReader.nrOfCol();++icol){
1069 // if(icol%down_opt[0])
1070  // continue;
1071 
1072 
1073  // if(rand_opt[0]>0){
1074  // gsl_rng* r=stat.getRandomGenerator(time(NULL));
1075  // //todo: init random number generator using time...
1076  // if(verbose_opt[0])
1077  // std::cout << "generating " << rand_opt[0] << " random numbers: " << std::endl;
1078  // for(unsigned int i=0;i<rand_opt[0];++i)
1079  // std::cout << i << " " << stat.getRandomValue(r,randdist_opt[0],randa_opt[0],randb_opt[0]) << std::endl;
1080  // }
1081 
1082  // imgreg.setDown(down_opt[0]);
1083  // imgreg.setThreshold(threshold_opt[0]);
1084  // double c0=0;//offset
1085  // double c1=1;//scale
1086  // double err=uncertNodata_opt[0];//start with high initial value in case we do not have first ob err=imgreg.getRMSE(imgReaderModel1,imgReader,c0,c1,verbose_opt[0]);
1087 
1088  // int nband=band_opt.size();
1089  // if(band_opt[0]<0)
1090  // nband=imgReader.nrOfBand();
1091  // for(int iband=0;iband<nband;++iband){
1092  // unsigned short band_opt[iband]=(band_opt[0]<0)? iband : band_opt[iband];
1093 
1094  // if(minmax_opt[0]||min_opt[0]||max_opt[0]){
1095  // assert(band_opt[iband]<imgReader.nrOfBand());
1096  // 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]))){
1097  // double uli,ulj,lri,lrj;
1098  // imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
1099  // imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
1100  // imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue);
1101  // }
1102  // else
1103  // imgReader.getMinMax(minValue,maxValue,band_opt[iband],true);
1104  // if(minmax_opt[0])
1105  // std::cout << "-min " << minValue << " -max " << maxValue << " ";
1106  // else{
1107  // if(min_opt[0])
1108  // std::cout << "-min " << minValue << " ";
1109  // if(max_opt[0])
1110  // std::cout << "-max " << maxValue << " ";
1111  // }
1112  // }
1113  // }
1114  // if(relative_opt[0])
1115  // hist_opt[0]=true;
1116  // if(hist_opt[0]){
1117  // assert(band_opt[0]<imgReader.nrOfBand());
1118  // unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
1119  // std::vector<unsigned long int> output;
1120  // minValue=0;
1121  // maxValue=0;
1122  // //todo: optimize such that getMinMax is only called once...
1123  // imgReader.getMinMax(minValue,maxValue,band_opt[0]);
1124 
1125  // if(src_min_opt.size())
1126  // minValue=src_min_opt[0];
1127  // if(src_max_opt.size())
1128  // maxValue=src_max_opt[0];
1129  // unsigned long int nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0]);
1130  // std::cout.precision(10);
1131  // for(int bin=0;bin<nbin;++bin){
1132  // double binValue=0;
1133  // if(nbin==maxValue-minValue+1)
1134  // binValue=minValue+bin;
1135  // else
1136  // binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin;
1137  // std::cout << binValue << " ";
1138  // if(relative_opt[0])
1139  // std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
1140  // else
1141  // std::cout << static_cast<double>(output[bin]) << std::endl;
1142  // }
1143  // }
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
void setOffset(double theOffset, int band=0)
Set offset for a specific band when writing the raster data values. The scaling and offset are applie...
Definition: ImgRasterGdal.h:84
bool isNoData(double value) const
Check if value is nodata in this dataset.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
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.
int nrOfRow(void) const
Get the number of rows of this dataset.
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked, irrespective of no data values.
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
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
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 readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
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...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void setScale(double theScale, int band=0)
Set scale for a specific band when writing the raster data values. The scaling and offset are applied...
Definition: ImgRasterGdal.h:75
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98