23#include "base/Optionpk.h"
24#include "algorithms/StatFactory.h"
25#include "algorithms/ImgRegression.h"
80int main(
int argc,
char *argv[])
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);
92 Optionpk<short> down_opt(
"down",
"down",
"Down sampling factor (for raster sample datasets only). Can be used to create grid points", 1);
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)");
102 Optionpk<bool> median_opt(
"median",
"median",
"calculate median",
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);
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);
142 doProcess=input_opt.retrieveOption(argc,argv);
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);
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);
178 catch(
string predefinedString){
179 std::cout << predefinedString << std::endl;
184 cout <<
"Usage: pkstat -i input" << endl;
186 std::cout <<
"short option -h shows basic options only, use long option --help to show all options" << std::endl;
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]);
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]);
204 double minValue=(src_min_opt.size())? src_min_opt[0] : 0;
205 double maxValue=(src_max_opt.size())? src_max_opt[0] : 0;
207 double medianValue=0;
210 const char* pszMessage;
211 void* pProgressArg=NULL;
212 GDALProgressFunc pfnProgress=GDALTermProgress;
218 std::vector<double> histogramOutput;
223 if(scale_opt.size()){
224 while(scale_opt.size()<input_opt.size())
225 scale_opt.push_back(scale_opt[0]);
227 if(offset_opt.size()){
228 while(offset_opt.size()<input_opt.size())
229 offset_opt.push_back(offset_opt[0]);
231 if(input_opt.empty()){
232 std::cerr <<
"No image dataset provided (use option -i). Use --help for help information";
235 for(
int ifile=0;ifile<input_opt.size();++ifile){
237 imgReader.
open(input_opt[ifile]);
239 catch(std::string errorstring){
240 std::cout << errorstring << std::endl;
245 std::cout <<
" --input " << input_opt[ifile] <<
" ";
247 for(
int inodata=0;inodata<nodata_opt.size();++inodata)
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){
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]);
264 if(stat_opt[0]||mean_opt[0]||median_opt[0]||var_opt[0]||stdev_opt[0]){
266 vector<double> readBuffer;
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);
274 std::cout <<
"--mean " << meanValue <<
" ";
276 std::cout <<
"--median " << medianValue <<
" ";
278 std::cout <<
"--stdDev " << sqrt(varValue) <<
" ";
280 std::cout <<
"--var " << varValue <<
" ";
282 std::cout <<
"-min " << minValue <<
" -max " << maxValue <<
" --mean " << meanValue <<
" --stdDev " << sqrt(varValue) <<
" ";
286 assert(band_opt[iband]<imgReader.
nrOfBand());
287 GDALProgressFunc pfnProgress;
289 GDALRasterBand* rasterBand;
291 rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
293 std::cout <<
"-min " << minValue <<
" -max " << maxValue <<
" --mean " << meanValue <<
" --stdDev " << stdDev <<
" ";
296 if(minmax_opt[0]||min_opt[0]||max_opt[0]){
297 assert(band_opt[iband]<imgReader.
nrOfBand());
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);
306 imgReader.
getMinMax(minValue,maxValue,band_opt[iband]);
309 std::cout <<
"-min " << minValue <<
" -max " << maxValue <<
" ";
312 std::cout <<
"-min " << minValue <<
" ";
314 std::cout <<
"-max " << maxValue <<
" ";
318 if(histogram_opt[0]){
319 assert(band_opt[0]<imgReader.
nrOfBand());
320 nbin=(nbin_opt.size())? nbin_opt[0]:0;
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]);
331 cout <<
"number of valid pixels in image: " << imgReader.
getNvalid(band_opt[0]) << endl;
333 nsample+=imgReader.
getHistogram(histogramOutput,minValue,maxValue,nbin,band_opt[0],kde_opt[0]);
336 if(ifile==input_opt.size()-1){
337 std::cout.precision(10);
338 for(
int bin=0;bin<nbin;++bin){
340 if(nbin==maxValue-minValue+1)
341 binValue=minValue+bin;
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;
348 std::cout << static_cast<double>(histogramOutput[bin]) << std::endl;
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()){
360 if(src_max_opt.size()){
364 nbin=(nbin_opt.size())? nbin_opt[0]:0;
366 std::cerr <<
"Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
368 imgReader.
getMinMax(minX,maxX,band_opt[0]);
370 imgReader.
getMinMax(minY,maxY,band_opt[1]);
372 minValue=(minX<minY)? minX:minY;
373 maxValue=(maxX>maxY)? maxX:maxY;
375 std::cout <<
"min and max values: " << minValue <<
", " << maxValue << std::endl;
376 nbin=maxValue-minValue+1;
382 assert(band_opt[0]<imgReader.
nrOfBand());
383 assert(band_opt[1]<imgReader.
nrOfBand());
384 GDALProgressFunc pfnProgress;
386 GDALRasterBand* rasterBand;
390 rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData);
392 rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData);
394 double estimatedSize=1.0*imgReader.
getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
396 estimatedSize*=random_opt[0]/100.0;
397 sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2);
402 std::cout <<
"calculating 2d kernel density estimate with sigma " << sigma <<
" for bands " << band_opt[0] <<
" and " << band_opt[1] << std::endl;
404 std::cout <<
"calculating 2d histogram for bands " << band_opt[0] <<
" and " << band_opt[1] << std::endl;
405 std::cout <<
"nbin: " << nbin << std::endl;
409 vector< vector<double> > output;
412 imgReader.
getMinMax(minX,maxX,band_opt[0]);
414 imgReader.
getMinMax(minY,maxY,band_opt[1]);
417 std::ostringstream s;
418 s<<
"Error: could not calculate distribution (minX>=maxX)";
422 std::ostringstream s;
423 s<<
"Error: could not calculate distribution (minY>=maxY)";
427 for(
int i=0;i<nbin;++i){
428 output[i].resize(nbin);
429 for(
int j=0;j<nbin;++j)
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){
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){
446 double p=
static_cast<double>(rand())/(RAND_MAX);
451 if(imgReader.
isNoData(inputX[icol]))
453 if(imgReader.
isNoData(inputY[icol]))
456 if(inputX[icol]>=maxX)
458 else if(inputX[icol]<=minX)
461 binX=
static_cast<int>(
static_cast<double>(inputX[icol]-minX)/(maxX-minX)*nbin);
462 if(inputY[icol]>=maxY)
464 else if(inputY[icol]<=minX)
467 binY=
static_cast<int>(
static_cast<double>(inputY[icol]-minY)/(maxY-minY)*nbin);
469 assert(binX<output.size());
471 assert(binY<output[binX].size());
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){
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;
487 ++output[binX][binY];
491 cout <<
"number of valid pixels: " << nvalid << endl;
493 for(
int binX=0;binX<nbin;++binX){
495 for(
int binY=0;binY<nbin;++binY){
497 if(nbin==maxX-minX+1)
500 binValueX=minX+
static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
502 if(nbin==maxY-minY+1)
505 binValueY=minY+
static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
507 double value=
static_cast<double>(output[binX][binY]);
512 cout << binValueX <<
" " << binValueY <<
" " << value << std::endl;
518 if(reg_opt[0]&&input_opt.size()<2){
519 if(band_opt.size()<2)
521 imgreg.setDown(down_opt[0]);
522 imgreg.setThreshold(random_opt[0]);
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;
528 if(regerr_opt[0]&&input_opt.size()<2){
529 if(band_opt.size()<2)
531 imgreg.setDown(down_opt[0]);
532 imgreg.setThreshold(random_opt[0]);
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;
538 if(rmse_opt[0]&&input_opt.size()<2){
539 if(band_opt.size()<2)
541 vector<double> xBuffer(imgReader.
nrOfCol());
542 vector<double> yBuffer(imgReader.
nrOfCol());
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];
558 if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[0]||yValue>src_max_opt[0])
561 double e=xValue-yValue;
567 double correctNorm=nValid;
570 std::cout <<
" -rmse " << sqrt(mse) << std::endl;
572 if(preg_opt[0]&&input_opt.size()<2){
573 if(band_opt.size()<2)
575 imgreg.setDown(down_opt[0]);
576 imgreg.setThreshold(random_opt[0]);
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;
651 if(reg_opt[0]&&(input_opt.size()>1)){
652 imgreg.setDown(down_opt[0]);
653 imgreg.setThreshold(random_opt[0]);
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]);
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]);
669 if(offset_opt.size())
670 imgReader1.setOffset(offset_opt[0],band_opt[0]);
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]);
678 for(
int inodata=0;inodata<nodata_opt.size();++inodata){
680 imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);
681 imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];
683 imgReader1.pushNoDataValue(nodata_opt[inodata]);
684 imgReader2.pushNoDataValue(nodata_opt[inodata]);
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;
692 if(preg_opt[0]&&(input_opt.size()>1)){
693 imgreg.setDown(down_opt[0]);
694 imgreg.setThreshold(random_opt[0]);
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]);
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]);
710 if(offset_opt.size())
711 imgReader1.setOffset(offset_opt[0],band_opt[0]);
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]);
719 for(
int inodata=0;inodata<nodata_opt.size();++inodata){
721 imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);
722 imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];
724 imgReader1.pushNoDataValue(nodata_opt[inodata]);
725 imgReader2.pushNoDataValue(nodata_opt[inodata]);
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;
733 if(regerr_opt[0]&&(input_opt.size()>1)){
734 imgreg.setDown(down_opt[0]);
735 imgreg.setThreshold(random_opt[0]);
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]);
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]);
751 if(offset_opt.size())
752 imgReader1.setOffset(offset_opt[0],band_opt[0]);
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]);
760 for(
int inodata=0;inodata<nodata_opt.size();++inodata){
762 imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);
763 imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];
765 imgReader1.pushNoDataValue(nodata_opt[inodata]);
766 imgReader2.pushNoDataValue(nodata_opt[inodata]);
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;
774 if(rmse_opt[0]&&(input_opt.size()>1)){
775 imgreg.setDown(down_opt[0]);
776 imgreg.setThreshold(random_opt[0]);
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]);
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]);
792 if(offset_opt.size())
793 imgReader1.setOffset(offset_opt[0],band_opt[0]);
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]);
801 for(
int inodata=0;inodata<nodata_opt.size();++inodata){
803 imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);
804 imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];
806 imgReader1.pushNoDataValue(nodata_opt[inodata]);
807 imgReader2.pushNoDataValue(nodata_opt[inodata]);
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;
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]);
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]);
829 if(offset_opt.size())
830 imgReader1.setOffset(offset_opt[0],band_opt[0]);
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]);
838 for(
int inodata=0;inodata<nodata_opt.size();++inodata){
840 imgReader1.GDALSetNoDataValue(nodata_opt[0],band_opt[0]);
841 imgReader2.GDALSetNoDataValue(nodata_opt[0]),band_opt[1];
843 imgReader1.pushNoDataValue(nodata_opt[inodata]);
844 imgReader2.pushNoDataValue(nodata_opt[inodata]);
847 imgReader1.getMinMax(minX,maxX,band_opt[0]);
848 imgReader2.getMinMax(minY,maxY,band_opt[1]);
851 cout <<
"minX: " << minX << endl;
852 cout <<
"maxX: " << maxX << endl;
853 cout <<
"minY: " << minY << endl;
854 cout <<
"maxY: " << maxY << endl;
857 if(src_min_opt.size()){
861 if(src_max_opt.size()){
866 nbin=(nbin_opt.size())? nbin_opt[0]:0;
868 std::cerr <<
"Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
872 imgReader1.getMinMax(minX,maxX,band_opt[0]);
874 imgReader2.getMinMax(minY,maxY,band_opt[1]);
876 minValue=(minX<minY)? minX:minY;
877 maxValue=(maxX>maxY)? maxX:maxY;
879 std::cout <<
"min and max values: " << minValue <<
", " << maxValue << std::endl;
880 nbin=maxValue-minValue+1;
886 GDALProgressFunc pfnProgress;
888 GDALRasterBand* rasterBand;
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);
897 double estimatedSize=1.0*imgReader.
getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
899 estimatedSize*=random_opt[0]/100.0;
900 sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2);
905 std::cout <<
"calculating 2d kernel density estimate with sigma " << sigma <<
" for datasets " << input_opt[0] <<
" and " << input_opt[1] << std::endl;
907 std::cout <<
"calculating 2d histogram for datasets " << input_opt[0] <<
" and " << input_opt[1] << std::endl;
908 std::cout <<
"nbin: " << nbin << std::endl;
911 vector< vector<double> > output;
914 imgReader1.getMinMax(minX,maxX,band_opt[0]);
916 imgReader2.getMinMax(minY,maxY,band_opt[1]);
919 std::ostringstream s;
920 s<<
"Error: could not calculate distribution (minX>=maxX)";
924 std::ostringstream s;
925 s<<
"Error: could not calculate distribution (minY>=maxY)";
929 cout <<
"minX: " << minX << endl;
930 cout <<
"maxX: " << maxX << endl;
931 cout <<
"minY: " << minY << endl;
932 cout <<
"maxY: " << maxY << endl;
935 for(
int i=0;i<nbin;++i){
936 output[i].resize(nbin);
937 for(
int j=0;j<nbin;++j)
942 vector<double> inputX(imgReader1.nrOfCol());
943 vector<double> inputY(imgReader2.nrOfCol());
951 for(
int irow=0;irow<imgReader1.nrOfRow();++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){
965 double p=
static_cast<double>(rand())/(RAND_MAX);
970 if(imgReader1.isNoData(inputX[icol]))
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]))
978 if(inputX[icol1]>=maxX)
980 else if(inputX[icol]<=minX)
983 binX=
static_cast<int>(
static_cast<double>(inputX[icol1]-minX)/(maxX-minX)*nbin);
984 if(inputY[icol2]>=maxY)
986 else if(inputY[icol2]<=minY)
989 binY=
static_cast<int>(
static_cast<double>(inputY[icol2]-minY)/(maxY-minY)*nbin);
991 assert(binX<output.size());
993 assert(binY<output[binX].size());
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){
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;
1010 ++output[binX][binY];
1016 cout <<
"number of valid pixels: " << nvalid << endl;
1017 for(
int binX=0;binX<nbin;++binX){
1019 for(
int binY=0;binY<nbin;++binY){
1021 if(nbin==maxX-minX+1)
1022 binValueX=minX+binX;
1024 binValueX=minX+
static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
1026 if(nbin==maxY-minY+1)
1027 binValueY=minY+binY;
1029 binValueY=minY+
static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
1030 double value=
static_cast<double>(output[binX][binY]);
1032 if(relative_opt[0]||kde_opt[0])
1033 value*=100.0/nvalid;
1035 cout << binValueX <<
" " << binValueY <<
" " << value << std::endl;
1044 if(!histogram_opt[0]||histogram2d_opt[0])
1045 std::cout << std::endl;
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...
int nrOfCol(void) const
Get the number of columns of this dataset.
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...
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...
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...