| 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
216 statfactory::StatFactory stat;
217 imgregression::ImgRegression imgreg;
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)
265 statfactory::StatFactory stat;
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// }
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::geo2image
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
Definition: ImgRasterGdal.cc:387
ImgRasterGdal::isNoData
bool isNoData(double value) const
Check if value is nodata in this dataset.
Definition: ImgRasterGdal.h:136
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::setScale
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
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgRasterGdal::setOffset
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
ImgRasterGdal::getRasterBand
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
Definition: ImgRasterGdal.cc:105
ImgRasterGdal::pushNoDataValue
int pushNoDataValue(double noDataValue)
Push a no data value for this dataset.
Definition: ImgRasterGdal.cc:492
ImgRasterGdal::GDALSetNoDataValue
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...
Definition: ImgRasterGdal.h:142
ImgRasterGdal::covers
bool covers(double x, double y) const
Check if a geolocation is covered by this dataset. Only the bounding box is checked,...
Definition: ImgRasterGdal.cc:443
ImgReaderGdal
Definition: ImgReaderGdal.h:37
ImgReaderGdal::getNvalid
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
Definition: ImgReaderGdal.cc:373
ImgReaderGdal::readData
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
ImgReaderGdal::getMinMax
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...
Definition: ImgReaderGdal.cc:161
ImgReaderGdal::close
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
Definition: ImgReaderGdal.cc:42
ImgReaderGdal::open
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
Definition: ImgReaderGdal.cc:36
ImgReaderGdal::readDataBlock
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...
Definition: ImgReaderGdal.h:220
ImgReaderGdal::getHistogram
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...
Definition: ImgReaderGdal.cc:257
Optionpk
Definition: Optionpk.h:107
imgregression::ImgRegression
Definition: ImgRegression.h:30
statfactory::StatFactory
Definition: StatFactory.h:43