pktools 2.6.7
Processing Kernel for geospatial data
pkstat.cc
1/**********************************************************************
2pkstat.cc: program to calculate basic statistics from raster dataset
3Copyright (C) 2008-2015 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 <iostream>
21#include <fstream>
22#include <math.h>
23#include "base/Optionpk.h"
24#include "algorithms/StatFactory.h"
25#include "algorithms/ImgRegression.h"
26/******************************************************************************/
78using namespace std;
79
80int 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// }
int nrOfRow(void) const
Get the number of rows of 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)
bool isNoData(double value) const
Check if value is nodata in this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
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
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
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
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...
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked,...
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
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
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 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...
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...