[pktools] 290/375: pkkalman with uncertainty band in observation input

Bas Couwenberg sebastic at xs4all.nl
Wed Dec 3 21:54:23 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 fb26c6fa17998d4dd9cc6fcdc4b822faf26448d1
Author: Pieter Kempeneers <kempenep at gmail.com>
Date:   Fri Jun 20 18:06:48 2014 +0200

    pkkalman with uncertainty band in observation input
---
 src/algorithms/ImgRegression.cc  |  2 +-
 src/apps/pkkalman.cc             | 64 +++++++++++++++++++++++++---------------
 src/imageclasses/ImgReaderGdal.h |  1 +
 src/imageclasses/ImgWriterGdal.h |  1 +
 4 files changed, 44 insertions(+), 24 deletions(-)

diff --git a/src/algorithms/ImgRegression.cc b/src/algorithms/ImgRegression.cc
index 8cb444d..b945206 100644
--- a/src/algorithms/ImgRegression.cc
+++ b/src/algorithms/ImgRegression.cc
@@ -78,6 +78,6 @@ double ImgRegression::getRMSE(const ImgReaderGdal& imgReader1, const ImgReaderGd
   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;
+    std::cout << "linear regression 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 0d53bcf..9d31cd9 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -50,8 +50,8 @@ int main(int argc,char **argv) {
   Optionpk<double> uncertModel_opt("um", "uncertmodel", "Multiply this value with std dev of first model image to obtain uncertainty of model",2);
   Optionpk<double> uncertObs_opt("uo", "uncertobs", "Uncertainty of valid observations",0);
   Optionpk<double> uncertNodata_opt("unodata", "uncertnodata", "Uncertainty in case of no-data values in observation", 10000);
-  Optionpk<double> regTime_opt("rt", "regtime", "Relative Weight for regression in time series", 0.5);
-  Optionpk<double> regSensor_opt("rs", "regsensor", "Relative Weight for regression in sensor series", 0.5);
+  // Optionpk<double> regTime_opt("rt", "regtime", "Relative Weight for regression in time series", 1.0);
+  // Optionpk<double> regSensor_opt("rs", "regsensor", "Relative Weight for regression in sensor series", 1.0);
   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.");
@@ -75,8 +75,8 @@ int main(int argc,char **argv) {
     uncertModel_opt.retrieveOption(argc,argv);
     uncertObs_opt.retrieveOption(argc,argv);
     uncertNodata_opt.retrieveOption(argc,argv);
-    regTime_opt.retrieveOption(argc,argv);
-    regSensor_opt.retrieveOption(argc,argv);
+    // regTime_opt.retrieveOption(argc,argv);
+    // regSensor_opt.retrieveOption(argc,argv);
     down_opt.retrieveOption(argc,argv);
     verbose_opt.retrieveOption(argc,argv);
   }
@@ -187,7 +187,6 @@ int main(int argc,char **argv) {
   double c1obs=0;
   double errObs=uncertNodata_opt[0];//start with high initial value in case we do not have first observation at time 0
 
-  //todo: map modindex to real time (e.g., Julian day)
   vector<int> relobsindex;
   // cout << tmodel_opt << endl;
   // cout << tobservation_opt << endl;
@@ -273,10 +272,14 @@ int main(int argc,char **argv) {
 	    double uncertObs=uncertObs_opt[0];
 	    if(uncertObsBuffer.size()>icol)
 	      uncertObs=uncertObsBuffer[icol];
-	    if(uncertObs>eps_opt[0]){
-	      double noemer=uncertObs*uncertObs+stdDev*stdDev;
-	      estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
-	      estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo:check! error?
+	    if(uncertModel_opt[0]*stdDev+uncertObs>eps_opt[0]){
+	      // double noemer=uncertObs*uncertObs+stdDev*stdDev;//todo: multiply stdDev with uncertModel_opt[0]
+	      // estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
+	      // estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo:check! error?
+	      double kalmanGain=uncertModel_opt[0]*stdDev/(uncertModel_opt[0]*stdDev+uncertObs);
+	      assert(kalmanGain<=1);
+	      estWriteBuffer[icol]=estReadBuffer[icol]+kalmanGain*(estWriteBuffer[icol]-estReadBuffer[icol]);
+	      uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev*(1-kalmanGain);
 	    }
 	    else{
 	      //no need to fill write buffer (already done in imgReaderObs.readData
@@ -329,7 +332,10 @@ int main(int argc,char **argv) {
       double c0mod=0;
       double c1mod=0;
 
-      double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+      if(verbose_opt[0])
+	  cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderModel2.getFileName() << endl;
+      double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod);
+      // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
 
       bool update=false;
       if(obsindex<relobsindex.size()){
@@ -343,6 +349,8 @@ int main(int argc,char **argv) {
 	imgReaderObs.getGeoTransform(geotransform);
 	imgReaderObs.setNoData(obsnodata_opt);
 	//calculate regression between model and observation
+	if(verbose_opt[0])
+	  cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderObs.getFileName() << endl;
 	errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
       }
       string input;
@@ -377,11 +385,9 @@ int main(int argc,char **argv) {
 	  double certNorm=(errMod*errMod+errObs*errObs);
 	  double certMod=errObs*errObs/certNorm;
 	  double certObs=errMod*errMod/certNorm;
-	  double regTime=(c0mod+c1mod*estValue)*certMod*regTime_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
-	  double regSensor=(c0obs+c1obs*estValue)*certObs*regSensor_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
+	  double regTime=(c0mod+c1mod*estValue)*certMod;
+	  double regSensor=(c0obs+c1obs*estValue)*certObs;
 	  estWriteBuffer[icol]=regTime+regSensor;
-	  // estWriteBuffer[icol]=(c0mod+c1mod*estValue)*certMod;
-	  // estWriteBuffer[icol]+=(c0obs+c1obs*estValue)*certObs;
 	  double totalUncertainty=0;
 	  if(errMod<eps_opt[0])
 	    totalUncertainty=errObs;
@@ -491,11 +497,20 @@ int main(int argc,char **argv) {
 	    double uncertObs=uncertObs_opt[0];
 	    if(uncertObsBuffer.size()>icol)
 	      uncertObs=uncertObsBuffer[icol];
-	    if(uncertObs>eps_opt[0]){
-	      double noemer=uncertObs*uncertObs+stdDev*stdDev;
-	      estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
-	      estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo: check error?
+	    if(uncertModel_opt[0]*stdDev+uncertObs>eps_opt[0]){
+	      // double noemer=uncertObs*uncertObs+stdDev*stdDev;//todo: multiply stdDev with uncertModel_opt[0]
+	      // estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
+	      // estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo:check! error?
+	      double kalmanGain=uncertModel_opt[0]*stdDev/(uncertModel_opt[0]*stdDev+uncertObs);
+	      assert(kalmanGain<=1);
+	      estWriteBuffer[icol]=estReadBuffer[icol]+kalmanGain*(estWriteBuffer[icol]-estReadBuffer[icol]);
+	      uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev*(1-kalmanGain);
 	    }
+	    // if(uncertObs>eps_opt[0]){
+	    //   double noemer=uncertObs*uncertObs+stdDev*stdDev;
+	    //   estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
+	    //   estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;//todo: check error?
+	    // }
 	    else{
 	      //no need to fill write buffer (already done in imgReaderObs.readData
 	      uncertWriteBuffer[icol]=uncertObs;
@@ -545,7 +560,10 @@ int main(int argc,char **argv) {
       double c0mod=0;
       double c1mod=0;
 
-      double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+      if(verbose_opt[0])
+	cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderModel2.getFileName() << endl;
+      double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod);
+      // double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
 
       bool update=false;
       if(obsindex<relobsindex.size()){
@@ -558,6 +576,8 @@ int main(int argc,char **argv) {
 	imgReaderObs.getGeoTransform(geotransform);
 	imgReaderObs.setNoData(obsnodata_opt);
 	//calculate regression between model and observation
+	if(verbose_opt[0])
+	  cout << "Calculating regression for " << imgReaderModel1.getFileName() << " " << imgReaderObs.getFileName() << endl;
 	errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
       }
       //prediction (also to fill cloudy pixels in update mode)
@@ -595,11 +615,9 @@ int main(int argc,char **argv) {
 	  double certNorm=(errMod*errMod+errObs*errObs);
 	  double certMod=errObs*errObs/certNorm;
 	  double certObs=errMod*errMod/certNorm;
-	  double regTime=(c0mod+c1mod*estValue)*certMod*regTime_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
-	  double regSensor=(c0obs+c1obs*estValue)*certObs*regSensor_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
+	  double regTime=(c0mod+c1mod*estValue)*certMod;
+	  double regSensor=(c0obs+c1obs*estValue)*certObs;
 	  estWriteBuffer[icol]=regTime+regSensor;
-	  // estWriteBuffer[icol]=(c0mod+c1mod*estValue)*certMod;
-	  // estWriteBuffer[icol]+=(c0obs+c1obs*estValue)*certObs;
 	  double totalUncertainty=0;
 	  if(errMod<eps_opt[0])
 	    totalUncertainty=errObs;
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index 5b68630..5d4fceb 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -38,6 +38,7 @@ public:
   ~ImgReaderGdal(void);
   void open(const std::string& filename);//, double magicX=1, double magicY=1);
   void close(void);
+  std::string getFileName() const {return m_filename;};
   int nrOfCol(void) const { return m_ncol;};
   int nrOfRow(void) const { return m_nrow;};
   int nrOfBand(void) const { return m_nband;};
diff --git a/src/imageclasses/ImgWriterGdal.h b/src/imageclasses/ImgWriterGdal.h
index a90f586..9cfc9a0 100644
--- a/src/imageclasses/ImgWriterGdal.h
+++ b/src/imageclasses/ImgWriterGdal.h
@@ -39,6 +39,7 @@ public:
   // void open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType="GTiff", const std::string& interleave="BAND", const std::string& compression="LZW", int magicX=1, int magicY=1);
   void open(const std::string& filename, int ncol, int nrow, int nband, const GDALDataType& dataType, const std::string& imageType, const std::vector<std::string>& options=std::vector<std::string>());
   void close(void);
+  std::string getFileName() const {return m_filename;};
   int nrOfCol(void) const { return m_ncol;};
   int nrOfRow(void) const { return m_nrow;};
   int nrOfBand(void) const { return m_nband;};

-- 
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