[pktools] 282/375: drafted, not tested pkkalman.cc
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:22 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit 0ae6a8ffdddd9a783a66361b3f44161ab8408f13
Author: user <user at osgeolive.(none)>
Date: Wed Jun 11 01:04:01 2014 +0200
drafted, not tested pkkalman.cc
---
src/algorithms/ImgRegression.cc | 5 +-
src/apps/pkkalman.cc | 186 ++++++++++++++++++++++++----------------
2 files changed, 117 insertions(+), 74 deletions(-)
diff --git a/src/algorithms/ImgRegression.cc b/src/algorithms/ImgRegression.cc
index 3ebf557..8cb444d 100644
--- a/src/algorithms/ImgRegression.cc
+++ b/src/algorithms/ImgRegression.cc
@@ -36,6 +36,7 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd
std::vector<double> buffer1;
std::vector<double> buffer2;
+ srand(time(NULL));
for(irow1=0;irow1<imgReader1.nrOfRow();++irow1){
if(irow1%m_down)
continue;
@@ -64,9 +65,8 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd
//check for nodata
double value1=rowBuffer1[icol1];
double value2=rowBuffer2[icol2];
- bool readValid=true;
if(imgReader1.isNoData(value1)||imgReader2.isNoData(value2))
- readValid=false;
+ continue;
buffer1.push_back(value1);
buffer2.push_back(value2);
@@ -76,6 +76,7 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd
}
statfactory::StatFactory stat;
double err=stat.linear_regression_err(buffer1,buffer2,c0,c1);
+
if(verbose)
std::cout << "linear regression model-model based on " << buffer1.size() << " points: " << c0 << "+" << c1 << " * x " << " with rmse: " << err << std::endl;
return err;
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index 38afbe5..6243185 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -34,11 +34,14 @@ int main(int argc,char **argv) {
Optionpk<string> model_opt("mod","model","model input datasets, e.g., MODIS (use: -mod model1 -mod model2 etc.");
Optionpk<string> observation_opt("obs","observation","observation input datasets, e.g., landsat (use: -obs obs1 -obs obs2 etc.");
Optionpk<int> tobservation_opt("tobs","tobservation","time sequence of observation input (sequence must have exact same length as observation input)");
- Optionpk<string> output_opt("o", "output", "Suffix for output image datasets");
+ Optionpk<string> output_opt("o", "output", "Output raster dataset");
Optionpk<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0).", 0);
Optionpk<double> modnodata_opt("modnodata", "modnodata", "invalid value for model input", 0);
Optionpk<double> obsnodata_opt("obsnodata", "obsnodata", "invalid value for observation input", 0);
- Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression", 90);
+ Optionpk<double> eps_opt("eps", "eps", "epsilon for non zero division", 0.00001);
+ Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression", 9);
+ Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image","GTiff",2);
+ Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0);
bool doProcess;//stop process when program was invoked with help option (-h --help)
@@ -50,6 +53,7 @@ int main(int argc,char **argv) {
threshold_opt.retrieveOption(argc,argv);
modnodata_opt.retrieveOption(argc,argv);
obsnodata_opt.retrieveOption(argc,argv);
+ eps_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
@@ -91,9 +95,25 @@ int main(int argc,char **argv) {
}
imgregression::ImgRegression imgreg;
- vector<ImgReaderGdal> imgReaderModel(model_opt.size());
- vector<ImgReaderGdal> imgReaderObs(observation_opt.size());
- vector<ImgWriterGdal> imgWriterPred(model_opt.size());
+ // vector<ImgReaderGdal> imgReaderModel(model_opt.size());
+ // vector<ImgReaderGdal> imgReaderObs(observation_opt.size());
+ ImgReaderGdal imgReaderModel1;
+ ImgReaderGdal imgReaderModel2;
+ ImgReaderGdal imgReaderObs;
+ ImgWriterGdal imgWriterEst;
+
+ imgReaderObs.open(observation_opt[0]);
+ int ncol=imgReaderObs.nrOfCol();
+ int nrow=imgReaderObs.nrOfRow();
+
+ string imageType=imgReaderObs.getImageType();
+ if(oformat_opt.size())//default
+ imageType=oformat_opt[0];
+ if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
+ string theInterleave="INTERLEAVE=";
+ theInterleave+=imgReaderObs.getInterleave();
+ option_opt.push_back(theInterleave);
+ }
int obsindex=0;
//forward step
@@ -102,86 +122,108 @@ int main(int argc,char **argv) {
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
double progress=0;
- srand(time(NULL));
imgreg.setDown(down_opt[0]);
imgreg.setThreshold(threshold_opt[0]);
+ double c0obs=0;
+ double c1obs=0;
+ double errObs=0;
for(int modindex=0;modindex<model_opt.size()-1;++modindex){
+ string output;
+ if(output_opt.size()==model_opt.size())
+ output=output_opt[modindex+1];
+ else{
+ ostringstream outputstream;
+ outputstream << output_opt[0] << "_" << modindex+1 << ".tif";
+ output=outputstream.str();
+ }
+
+ //two band output band0=estimation, band1=uncertainty
+ imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+
//calculate regression between two subsequence model inputs
- imgReaderModel[modindex].open(model_opt[modindex]);
- imgReaderModel[modindex+1].open(model_opt[modindex+1]);
- imgReaderModel[modindex].setNoData(modnodata_opt);
- imgReaderModel[modindex+1].setNoData(modnodata_opt);
+ imgReaderModel1.open(model_opt[modindex]);
+ imgReaderModel2.open(model_opt[modindex+1]);
//calculate regression
//we could re-use the points from second image from last run, but
//to keep it general, we must redo it (overlap might have changed)
+
pfnProgress(progress,pszMessage,pProgressArg);
- int icol1=0,irow1=0;
- vector<double> rowBuffer1(imgReaderModel[modindex].nrOfCol());
- vector<double> rowBuffer2(imgReaderModel[modindex+1].nrOfCol());
- vector<double> buffer1;
- vector<double> buffer2;
-
- double c0=0;
- double c1=0;
-
- // for(irow1=0;irow1<imgReaderModel[modindex].nrOfRow();++irow1){
- // if(irow1%down_opt[0])
- // continue;
- // icol1=0;
- // double icol2=0,irow2=0;
- // double geox=0,geoy=0;
- // imgReaderModel[modindex].readData(rowBuffer1,GDT_Float64,irow1);
- // imgReaderModel[modindex].image2geo(icol1,irow1,geox,geoy);
- // imgReaderModel[modindex+1].geo2image(geox,geoy,icol2,irow2);
- // icol2=static_cast<int>(icol2);
- // irow2=static_cast<int>(irow2);
- // imgReaderModel[modindex+1].readData(rowBuffer2,GDT_Float64,irow2);
- // for(icol1=0;icol1<imgReaderModel[modindex].nrOfCol();++icol1){
- // if(icol1%down_opt[0])
- // continue;
- // if(threshold_opt[0]>0){//percentual value
- // double p=static_cast<double>(rand())/(RAND_MAX);
- // p*=100.0;
- // if(p>threshold_opt[0])
- // continue;//do not select for now, go to next column
- // }
- // imgReaderModel[modindex].image2geo(icol1,irow1,geox,geoy);
- // imgReaderModel[modindex+1].geo2image(geox,geoy,icol2,irow2);
- // icol2=static_cast<int>(icol2);
- // irow2=static_cast<int>(irow2);
- // //check for nodata
- // double valmod1=rowBuffer1[icol1];
- // double valmod2=rowBuffer2[icol2];
- // bool readValid=true;
- // for(int vband=0;vband<bndmodnodata_opt.size();++vband){
- // if(modnodata_opt.size()>vband){
- // if(valmod1==modnodata_opt[vband] || valmod2==modnodata_opt[vband]){
- // readValid=false;
- // break;
- // }
- // }
- // }
- // buffer1.push_back(valmod1);
- // buffer2.push_back(valmod2);
- // if(verbose_opt[0]>1)
- // cout << geox << " " << geoy << " " << icol1 << " " << irow1 << " " << icol2 << " " << irow2 << " " << buffer1.back() << " " << buffer2.back() << endl;
- // }
- // }
-
- double err=imgreg.getRMSE(imgReaderModel[modindex],imgReaderModel[modindex+1],c0,c1,verbose_opt[0]);
-
- if(tobservation_opt[obsindex]==modindex){//update
- imgReaderObs[obsindex].open(observation_opt[obsindex]);
- //calculate regression
- //set ImgWriterPred[modindex]=imgReaderObs[obindex]
- //calculate regression between model and observation
+ double c0mod=0;
+ double c1mod=0;
+
+ double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+ bool update=(tobservation_opt[obsindex]==modindex);
+ if(update){//update
+ imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.setNoData(obsnodata_opt);
+ //calculate regression between model and observation
+ errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ }
+ //prediction (also to fill cloudy pixels in update mode)
+ string input;
+ if(output_opt.size()==model_opt.size())
+ //todo: initialize output_opt[0] with model[0]
+ input=output_opt[modindex];
+ else{
+ ostringstream outputstring;
+ outputstring << output_opt[0] << "_" << modindex << ".tif";
+ input=outputstring.str();
+ }
+ ImgReaderGdal imgReaderEst(input);
+
+ vector<double> obsBuffer;
+ vector<double> estReadBuffer;
+ vector<double> estWriteBuffer;
+ vector<double> uncertWriteBuffer;
+ for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
+ assert(irow<imgReaderEst.nrOfRow());
+ imgReaderEst.readData(estReadBuffer,GDT_Float64,irow);
+ if(update){
+ imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow);
+ //todo: write uncertainty image: 0 if observation,
+ }
+ for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
+ if(update){//check for nodata in observation
+ if(!imgReaderObs.isNoData(estWriteBuffer[icol])){
+ uncertWriteBuffer[icol]=0;
+ continue;//no need to estimate
+ }
+ }
+ //predict
+ double estValue=estReadBuffer[icol];
+ double certNorm=(errMod*errMod+errObs*errObs);
+ double certMod=errObs*errObs/certNorm;
+ double certObs=errMod*errMod/certNorm;
+ estWriteBuffer[icol]=(c0mod+c1mod*estValue)*certMod;
+ estWriteBuffer[icol]+=(c0obs+c1obs*estValue)*certObs;
+
+ double totalUncertainty=0;
+ if(errMod<eps_opt[0])
+ totalUncertainty=errObs;
+ else if(errObs<eps_opt[0])
+ totalUncertainty=errMod;
+ else{
+ totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
+ totalUncertainty=sqrt(1.0/totalUncertainty);
+ }
+ uncertWriteBuffer[icol]=totalUncertainty;
+ }
+ imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
+ imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
}
- else{//prediction
+
+ imgWriterEst.close();
+ imgReaderEst.close();
+
+ if(update){
+ imgReaderObs.close();
+ ++obsindex;
}
- ++obsindex;
+ imgReaderModel1.close();
+ imgReaderModel2.close();
}
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list