[pktools] 286/375: redesigned forward kalman filter, todo: backward

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 9d9621cdf4906993e9346c8c269a53081df15cef
Author: user <user at osgeolive.(none)>
Date:   Fri Jun 13 01:43:59 2014 +0200

    redesigned forward kalman filter, todo: backward
---
 src/apps/pkkalman.cc | 231 +++++++++++++++++++++++++++++----------------------
 1 file changed, 133 insertions(+), 98 deletions(-)

diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index 886d7a1..73959d0 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -27,6 +27,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include "algorithms/StatFactory.h"
 #include "algorithms/ImgRegression.h"
 
+				    //todo: keep original resolution of coarse model raster dataset
 using namespace std;
 /*------------------
   Main procedure
@@ -45,7 +46,10 @@ int main(int argc,char **argv) {
   Optionpk<double> obsnodata_opt("obsnodata", "obsnodata", "invalid value for observation input", 0);
   Optionpk<double> eps_opt("eps", "eps", "epsilon for non zero division", 0.00001);
   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<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.");
@@ -66,7 +70,10 @@ int main(int argc,char **argv) {
     obsnodata_opt.retrieveOption(argc,argv);
     eps_opt.retrieveOption(argc,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);
     down_opt.retrieveOption(argc,argv);
     verbose_opt.retrieveOption(argc,argv);
   }
@@ -180,7 +187,7 @@ int main(int argc,char **argv) {
   for(int tindex=0;tindex<tobservation_opt.size();++tindex){
     vector<int>::iterator modit;
     modit=lower_bound(tmodel_opt.begin(),tmodel_opt.end(),tobservation_opt[tindex]);
-    int relpos=modit-tmodel_opt.begin()-2;
+    int relpos=modit-tmodel_opt.begin()-1;
     if(relpos<0)
       relpos=0;
     relobsindex.push_back(relpos);
@@ -208,57 +215,74 @@ int main(int argc,char **argv) {
     imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
 
 
+    imgReaderModel1.open(model_opt[0]);
+    //calculate standard deviation of image to serve as model uncertainty
+    GDALRasterBand* rasterBand;
+    rasterBand=imgReaderModel1.getRasterBand(0);
+    double minValue, maxValue, meanValue, stdDev;
+    void* pProgressData;
+    rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
     if(relobsindex[0]>0){//initialize output_opt[0] as model[0]
       //write first model as output
-      imgReaderModel1.open(model_opt[0]);
-      //calculate standard deviation of image to serve as model uncertainty
-      GDALRasterBand* rasterBand;
-      rasterBand=imgReaderModel1.getRasterBand(0);
-      double minValue, maxValue, meanValue, stdDev;
-      void* pProgressData;
-      rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
-
-      for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
-	vector<double> buffer;
-	imgReaderModel1.readData(buffer,GDT_Float64,irow);
-	imgWriterEst.writeData(buffer,GDT_Float64,irow,0);
-	for(int icol=0;icol<imgWriterEst.nrOfCol();++icol)
-	  buffer[icol]=uncertModel_opt[0]*stdDev;
-	imgWriterEst.writeData(buffer,GDT_Float64,irow,1);
+      for(int irow=0;irow<nrow;++irow){
+	vector<double> estReadBuffer;
+	vector<double> uncertWriteBuffer(ncol);
+	imgReaderModel1.readData(estReadBuffer,GDT_Float64,irow);
+	imgWriterEst.writeData(estReadBuffer,GDT_Float64,irow,0);
+	for(int icol=0;icol<ncol;++icol)
+	  uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
+	imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
       }
-      imgReaderModel1.close();
-      imgWriterEst.close();
     }
     else{//we have an observation at time 0
       imgReaderObs.open(observation_opt[0]);
       imgReaderObs.setNoData(obsnodata_opt);
-      for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
-	vector<double> buffer;
-	imgReaderObs.readData(buffer,GDT_Float64,irow);
-	imgWriterEst.writeData(buffer,GDT_Float64,irow,0);
-	for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
-	  if(imgReaderObs.isNoData(buffer[icol]))
-	    buffer[icol]=uncertNodata_opt[0];
-	  else
-	    buffer[icol]=0;
+      for(int irow=0;irow<nrow;++irow){
+	vector<double> estReadBuffer;
+	imgReaderModel1.readData(estReadBuffer,GDT_Float64,irow);
+	vector<double> estWriteBuffer(ncol);
+	vector<double> uncertWriteBuffer(ncol);
+	imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
+	for(int icol=0;icol<ncol;++icol){
+	  if(imgReaderObs.isNoData(estWriteBuffer[icol])){
+	    estWriteBuffer[icol]=estReadBuffer[icol];
+	    uncertWriteBuffer[icol]=uncertModel_opt[0]*stdDev;
+	  }
+	  else{
+	    if(uncertObs_opt[0]>eps_opt[0]){
+	      double noemer=uncertObs_opt[0]*uncertObs_opt[0]+stdDev*stdDev;
+	      estWriteBuffer[icol]*=uncertModel_opt[0]*stdDev*stdDev/noemer;
+	      estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs_opt[0]*uncertObs_opt[0]/noemer;
+	    }
+	    else{
+	      //no need to fill write buffer (already done in imgReaderObs.readData
+	      uncertWriteBuffer[icol]=uncertObs_opt[0];
+	    }
+	  }
 	}
-	imgWriterEst.writeData(buffer,GDT_Float64,irow,1);
+	imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
+	imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
       }
       imgReaderObs.close();
-      imgWriterEst.close();
+      ++obsindex;
     }
+    imgReaderModel1.close();
+    imgWriterEst.close();
 
-    for(int modindex=0;modindex<model_opt.size()-1;++modindex){
+    for(int modindex=1;modindex<model_opt.size();++modindex){
       if(verbose_opt[0]){
-	cout << "processing time " << modindex << endl;
-	cout << "last observation " << relobsindex[obsindex] << endl;
+	cout << "processing time " << tmodel_opt[modindex] << endl;
+	if(obsindex<relobsindex.size())
+	  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+	else
+	  cout << "There is no next observation" << endl;
       }
       string output;
       if(outputfw_opt.size()==model_opt.size())
-	output=outputfw_opt[modindex+1];
+	output=outputfw_opt[modindex];
       else{
 	ostringstream outputstream;
-	outputstream << outputfw_opt[0] << "_" << tmodel_opt[modindex+1] << ".tif";
+	outputstream << outputfw_opt[0] << "_" << tmodel_opt[modindex] << ".tif";
 	// outputstream << output_opt[0] << "_" << modindex+1 << ".tif";
 	output=outputstream.str();
       }
@@ -268,8 +292,8 @@ int main(int argc,char **argv) {
       imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
 
       //calculate regression between two subsequence model inputs
-      imgReaderModel1.open(model_opt[modindex]);
-      imgReaderModel2.open(model_opt[modindex+1]);
+      imgReaderModel1.open(model_opt[modindex-1]);
+      imgReaderModel2.open(model_opt[modindex]);
       //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)
@@ -280,22 +304,25 @@ int main(int argc,char **argv) {
 
       double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
 
-      bool update=(relobsindex[obsindex]==modindex);
+      bool update=false;
+      if(obsindex<relobsindex.size()){
+	update=(relobsindex[obsindex]==modindex);
+      }
       if(update){
 	if(verbose_opt[0])
-	  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " ***" << endl;
+	  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
+
 	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(outputfw_opt.size()==model_opt.size())
-	input=outputfw_opt[modindex];
+	input=outputfw_opt[modindex-1];
       else{
 	ostringstream outputstream;
-	outputstream << outputfw_opt[0] << "_" << tmodel_opt[modindex] << ".tif";
+	outputstream << outputfw_opt[0] << "_" << tmodel_opt[modindex-1] << ".tif";
 	input=outputstream.str();
       }
       ImgReaderGdal imgReaderEst(input);
@@ -310,42 +337,34 @@ int main(int argc,char **argv) {
 	imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
 	imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
 	if(update){
-	  imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow);
+	  imgReaderObs.readData(obsBuffer,GDT_Float64,irow);
 	}
 	for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
-	  double kalmanGain=0;
 	  double estValue=estReadBuffer[icol];
 	  double uncertValue=uncertReadBuffer[icol];
-	  if(update){//check for nodata in observation
-	    if(imgReaderObs.isNoData(estWriteBuffer[icol])){
-	      kalmanGain=0;
-	      // uncertWriteBuffer[icol]=0;
-	    }
-	    else
-	      kalmanGain=1;
-	    //todo: introduce gains between 0 and 1 in case of uncertain observations (SLC, potentially cloudy, etc.
-	    estWriteBuffer[icol]*=kalmanGain;
-	    estWriteBuffer[icol]+=(1-kalmanGain)*estReadBuffer[icol];
-	    uncertWriteBuffer[icol]=uncertReadBuffer[icol]*(1-kalmanGain);
-	  }
-	  else{//time update
-	    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;
+	  //time update
+	  double certNorm=(errMod*errMod+errObs*errObs);
+	  double certMod=errObs*errObs/certNorm;
+	  double certObs=errMod*errMod/certNorm;
+	  estWriteBuffer[icol]=(c0mod+c1mod*estValue)*certMod*regTime_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
+	  estWriteBuffer[icol]+=(c0obs+c1obs*estValue)*certObs*regSensor_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
 	  
-	    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+uncertReadBuffer[icol];
+	  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+uncertReadBuffer[icol];
+	  if(update){
+	    double kalmanGain=1;
+	    if((uncertWriteBuffer[icol]+uncertObs_opt[0])>eps_opt[0])
+	      kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs_opt[0]);
+	    estWriteBuffer[icol]+=kalmanGain*(obsBuffer[icol]-estWriteBuffer[icol]);
+	    uncertWriteBuffer[icol]*=(1-kalmanGain);
 	  }
 	}
 	imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
@@ -417,7 +436,7 @@ int main(int argc,char **argv) {
 	  if(imgReaderObs.isNoData(buffer[icol]))
 	    buffer[icol]=uncertNodata_opt[0];
 	  else
-	    buffer[icol]=0;
+	    buffer[icol]=uncertObs_opt[0];
 	}
 	imgWriterEst.writeData(buffer,GDT_Float64,irow,1);
       }
@@ -427,8 +446,11 @@ int main(int argc,char **argv) {
 
     for(int modindex=model_opt.size()-1;modindex>0;--modindex){
       if(verbose_opt[0]){
-	cout << "processing time " << modindex << endl;
-	cout << "last observation " << relobsindex[obsindex] << endl;
+	cout << "processing time " << tmodel_opt[modindex] << endl;
+	if(obsindex<relobsindex.size())
+	  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+	else
+	  cout << "There is no next observation" << endl;
       }
       string output;
       if(outputbw_opt.size()==model_opt.size())
@@ -457,10 +479,13 @@ int main(int argc,char **argv) {
 
       double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
 
-      bool update=(relobsindex[obsindex]==modindex);
+      bool update=false;
+      if(obsindex<relobsindex.size()){
+	update=(relobsindex[obsindex]==modindex);
+      }
       if(update){
 	if(verbose_opt[0])
-	  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " ***" << endl;
+	  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
 	imgReaderObs.open(observation_opt[obsindex]);
 	imgReaderObs.setNoData(obsnodata_opt);
 	//calculate regression between model and observation
@@ -498,9 +523,10 @@ int main(int argc,char **argv) {
 	      kalmanGain=0;
 	      // uncertWriteBuffer[icol]=0;
 	    }
-	    else
+	    else{
 	      kalmanGain=1;
 	    //todo: introduce gains between 0 and 1 in case of uncertain observations (SLC, potentially cloudy, etc.
+	    }
 	    estWriteBuffer[icol]*=kalmanGain;
 	    estWriteBuffer[icol]+=(1-kalmanGain)*estReadBuffer[icol];
 	    uncertWriteBuffer[icol]=uncertReadBuffer[icol]*(1-kalmanGain);
@@ -510,8 +536,8 @@ int main(int argc,char **argv) {
 	    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;
+	    estWriteBuffer[icol]=(c0mod+c1mod*estValue)*certMod*regTime_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
+	    estWriteBuffer[icol]+=(c0obs+c1obs*estValue)*certObs*regSensor_opt[0]/(regTime_opt[0]+regSensor_opt[0]);
 	  
 	    double totalUncertainty=0;
 	    if(errMod<eps_opt[0])
@@ -549,7 +575,11 @@ int main(int argc,char **argv) {
       cout << "Running smooth model" << endl;
     for(int modindex=0;modindex<model_opt.size();++modindex){
       if(verbose_opt[0]){
-	cout << "processing time " << modindex << endl;
+	cout << "processing time " << tmodel_opt[modindex] << endl;
+	if(obsindex<relobsindex.size())
+	  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+	else
+	  cout << "There is no next observation" << endl;
       }
       string output;
       if(outputfb_opt.size()==model_opt.size())
@@ -596,11 +626,14 @@ int main(int argc,char **argv) {
       vector<double> estWriteBuffer(ncol);
       vector<double> uncertWriteBuffer(ncol);
 
-      bool update=(relobsindex[obsindex]==modindex);
+      bool update=false;
+      if(obsindex<relobsindex.size()){
+	update=(relobsindex[obsindex]==modindex);
+      }
 
       if(update){
 	if(verbose_opt[0])
-	  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " ***" << endl;
+	  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
 	imgReaderObs.open(observation_opt[obsindex]);
 	imgReaderObs.setNoData(obsnodata_opt);
 	//calculate regression between model and observation
@@ -635,28 +668,30 @@ int main(int argc,char **argv) {
 	    else if(imgReaderObs.isNoData(estWriteBuffer[icol]))
 	      uncertObs=uncertNodata_opt[0];
 	    else
-	      uncertObs=0;
+	      uncertObs=uncertObs_opt[0];
 	  }
 
 	  double noemer=(C+D);
 	  //todo: consistently check for division by zero...
-	  if(noemer<eps_opt[0])
-	    estWriteBuffer[icol]=obsnodata_opt[0];
+	  if(noemer<eps_opt[0]){//simple average if both uncertainties are ~>0
+	    estWriteBuffer[icol]=0.5*(A+B);
+	    uncertWriteBuffer[icol]=uncertObs_opt[0];
+	  }
 	  else{
 	    estWriteBuffer[icol]=(A*D+B*C)/noemer;
+	    double P=0;
+	    if(C>eps_opt[0])
+	      P+=1.0/C;
+	    if(D>eps_opt[0])
+	      P+=1.0/D;
+	    if(uncertObs*uncertObs>eps_opt[0])
+	      P-=1.0/uncertObs/uncertObs;
+	    if(P>eps_opt[0])
+	      P=sqrt(1.0/P);
+	    else
+	      P=0;
+	    uncertWriteBuffer[icol]=P;
 	  }
-	  double P=0;
-	  if(C>eps_opt[0])
-	    P+=1.0/C;
-	  if(D>eps_opt[0])
-	    P+=1.0/D;
-	  if(uncertObs*uncertObs>eps_opt[0])
-	    P-=1.0/uncertObs/uncertObs;
-	  if(P>eps_opt[0])
-	    P=sqrt(1.0/P);
-	  else
-	    P=0;
-	  uncertWriteBuffer[icol]=P;
 	}
 	imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
 	imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);

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