[pktools] 284/375: implemented both forward and backward kalman filter
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 8f2e83e32dbfc11c91e698a2d5cc3e60e583d9eb
Author: user <user at osgeolive.(none)>
Date: Thu Jun 12 00:13:13 2014 +0200
implemented both forward and backward kalman filter
---
src/apps/pkkalman.cc | 495 ++++++++++++++++++++++++++++++++++++---------------
1 file changed, 356 insertions(+), 139 deletions(-)
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index 75ada00..f858826 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -19,6 +19,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include <sstream>
#include <vector>
+#include <algorithm>
#include "base/Optionpk.h"
#include "base/Vector2d.h"
#include "imageclasses/ImgReaderGdal.h"
@@ -31,15 +32,17 @@ using namespace std;
Main procedure
----------------*/
int main(int argc,char **argv) {
+ Optionpk<string> direction_opt("dir","direction","direction to run model (forward|backward)","forward");
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<int> tmodel_opt("tmod","tmodel","time sequence of model input. Sequence must have exact same length as model input. Leave empty to have default sequence 0,1,2,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", "Output raster dataset");
Optionpk<float> threshold_opt("th", "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<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");
+ Optionpk<double> uncertModel_opt("um", "uncertmodel", "Multiply this value with std dev of first model image to obtain uncertainty of model",2);
Optionpk<double> uncertNodata_opt("unodata", "uncertnodata", "Uncertainty in case of no-data values in observation", 10000);
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);
@@ -48,8 +51,10 @@ int main(int argc,char **argv) {
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
- doProcess=model_opt.retrieveOption(argc,argv);
+ doProcess=direction_opt.retrieveOption(argc,argv);
+ model_opt.retrieveOption(argc,argv);
observation_opt.retrieveOption(argc,argv);
+ tmodel_opt.retrieveOption(argc,argv);
tobservation_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
@@ -84,6 +89,15 @@ int main(int argc,char **argv) {
errorStream << "Error: sequence of models should be larger than observations" << endl;
throw(errorStream.str());
}
+ if(tmodel_opt.size()!=model_opt.size()){
+ if(tmodel_opt.empty())
+ cout << "Warning: time sequence is not provided, self generating time sequence from 0 to " << model_opt.size() << endl;
+ else
+ cout << "Warning: time sequence provided (" << tmodel_opt.size() << ") does not match number of model raster datasets (" << model_opt.size() << ")" << endl;
+ tmodel_opt.clear();
+ for(int tindex=0;tindex<model_opt.size();++tindex)
+ tmodel_opt.push_back(tindex);
+ }
if(tobservation_opt.size()!=observation_opt.size()){
errorStream << "Error: time sequence for observation must match size of observation dataset" << endl;
throw(errorStream.str());
@@ -120,7 +134,6 @@ int main(int argc,char **argv) {
}
int obsindex=0;
- //forward step
const char* pszMessage;
void* pProgressArg=NULL;
@@ -134,171 +147,375 @@ 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
- //initialization
- string output;
- if(output_opt.size()==model_opt.size())
- output=output_opt[0];
- else{
- ostringstream outputstream;
- outputstream << output_opt[0] << "_0.tif";
- output=outputstream.str();
- }
- imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
- imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
- if(tobservation_opt[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);
- }
- 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;
- }
- imgWriterEst.writeData(buffer,GDT_Float64,irow,1);
- }
- imgReaderObs.close();
- imgWriterEst.close();
- }
-
//todo: map modindex to real time (e.g., Julian day)
+ vector<int> relobsindex;
+ // cout << tmodel_opt << endl;
+ // cout << tobservation_opt << endl;
- for(int modindex=0;modindex<model_opt.size()-1;++modindex){
+ 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;
+ if(relpos<0)
+ relpos=0;
+ relobsindex.push_back(relpos);
if(verbose_opt[0])
- cout << "processing time " << modindex << endl;
+ cout << "tobservation_opt[tindex] " << tobservation_opt[tindex] << " " << relobsindex.back() << endl;
+ }
+
+ if(direction_opt[0]=="forward"){
+ ///////////////////////////// forward model /////////////////////////
+ obsindex=0;
+ if(verbose_opt[0])
+ cout << "Running forward model" << endl;
+ //initialization
string output;
if(output_opt.size()==model_opt.size())
- output=output_opt[modindex+1];
+ output=output_opt[0];
else{
ostringstream outputstream;
- outputstream << output_opt[0] << "_" << modindex+1 << ".tif";
+ outputstream << output_opt[0] << "_" << tmodel_opt[0] << ".tif";
output=outputstream.str();
}
-
- //two band output band0=estimation, band1=uncertainty
+ if(verbose_opt[0])
+ cout << "Opening image " << output << " for writing " << endl;
imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
- //calculate regression between two subsequence model inputs
- 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);
- double c0mod=0;
- double c1mod=0;
- double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+ 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);
- bool update=(tobservation_opt[obsindex]==modindex);
- if(update){
- if(verbose_opt[0])
- cout << "***update " << tobservation_opt[obsindex] << " = " << modindex << " ***" << endl;
- imgReaderObs.open(observation_opt[obsindex]);
+ 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);
+ }
+ imgReaderModel1.close();
+ imgWriterEst.close();
+ }
+ else{//we have an observation at time 0
+ imgReaderObs.open(observation_opt[0]);
imgReaderObs.setNoData(obsnodata_opt);
- //calculate regression between model and observation
- errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
+ 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;
+ }
+ imgWriterEst.writeData(buffer,GDT_Float64,irow,1);
+ }
+ imgReaderObs.close();
+ imgWriterEst.close();
+ }
+
+ for(int modindex=0;modindex<model_opt.size()-1;++modindex){
+ if(verbose_opt[0]){
+ cout << "processing time " << modindex << endl;
+ cout << "last observation " << relobsindex[obsindex] << endl;
+ }
+ string output;
+ if(output_opt.size()==model_opt.size())
+ output=output_opt[modindex+1];
+ else{
+ ostringstream outputstream;
+ outputstream << output_opt[0] << "_" << tmodel_opt[modindex+1] << ".tif";
+ // 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);
+ imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
+
+ //calculate regression between two subsequence model inputs
+ 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);
+ double c0mod=0;
+ double c1mod=0;
+
+ double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+
+ bool update=(relobsindex[obsindex]==modindex);
+ if(update){
+ if(verbose_opt[0])
+ cout << "***update " << relobsindex[obsindex] << " = " << modindex << " ***" << 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(output_opt.size()==model_opt.size())
+ input=output_opt[modindex];
+ else{
+ ostringstream outputstream;
+ outputstream << output_opt[0] << "_" << tmodel_opt[modindex] << ".tif";
+ input=outputstream.str();
+ }
+ ImgReaderGdal imgReaderEst(input);
+
+ vector<double> obsBuffer;
+ vector<double> estReadBuffer;
+ vector<double> uncertReadBuffer;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
+ assert(irow<imgReaderEst.nrOfRow());
+ imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
+ imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
+ if(update){
+ imgReaderObs.readData(estWriteBuffer,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;
+
+ 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];
+ }
+ }
+ imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
+ imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
+ progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+
+ imgWriterEst.close();
+ imgReaderEst.close();
+
+ if(update){
+ imgReaderObs.close();
+ ++obsindex;
+ }
+ imgReaderModel1.close();
+ imgReaderModel2.close();
}
- //prediction (also to fill cloudy pixels in update mode)
- string input;
+ }
+ else{
+ obsindex=relobsindex.size()-1;
+ ///////////////////////////// backward model /////////////////////////
+ if(verbose_opt[0])
+ cout << "Running backward model" << endl;
+ //initialization
+ string output;
if(output_opt.size()==model_opt.size())
- input=output_opt[modindex];
+ output=output_opt.back();
else{
- ostringstream outputstring;
- outputstring << output_opt[0] << "_" << modindex << ".tif";
- input=outputstring.str();
+ ostringstream outputstream;
+ outputstream << output_opt[0] << "_" << tmodel_opt.back() << ".tif";
+ output=outputstream.str();
}
- ImgReaderGdal imgReaderEst(input);
-
- vector<double> obsBuffer;
- vector<double> estReadBuffer;
- vector<double> uncertReadBuffer;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
- assert(irow<imgReaderEst.nrOfRow());
- imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
- imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
- if(update){
- imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow);
+ if(verbose_opt[0])
+ cout << "Opening image " << output << " for writing " << endl;
+ imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+ imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
+
+ //todo: check if correct...
+ if(relobsindex.back()<model_opt.size()-1){//initialize output_opt[0] as model[0]
+ //write last model as output
+ imgReaderModel1.open(model_opt.back());
+ //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 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;
- }
+ imgReaderModel1.close();
+ imgWriterEst.close();
+ }
+ else{//we have an observation at end time
+ imgReaderObs.open(observation_opt.back());
+ 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
- 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);
+ buffer[icol]=0;
+ }
+ imgWriterEst.writeData(buffer,GDT_Float64,irow,1);
+ }
+ imgReaderObs.close();
+ imgWriterEst.close();
+ }
+
+ 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;
+ }
+ string output;
+ if(output_opt.size()==model_opt.size())
+ output=output_opt[modindex-1];
+ else{
+ ostringstream outputstream;
+ outputstream << output_opt[0] << "_" << tmodel_opt[modindex-1] << ".tif";
+ // 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);
+ imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
+
+ //calculate regression between two subsequence model inputs
+ 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);
+ double c0mod=0;
+ double c1mod=0;
+
+ double errMod=imgreg.getRMSE(imgReaderModel1,imgReaderModel2,c0mod,c1mod,verbose_opt[0]);
+
+ bool update=(relobsindex[obsindex]==modindex);
+ if(update){
+ if(verbose_opt[0])
+ cout << "***update " << relobsindex[obsindex] << " = " << modindex << " ***" << 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(output_opt.size()==model_opt.size())
+ input=output_opt[modindex];
+ else{
+ ostringstream outputstream;
+ outputstream << output_opt[0] << "_" << tmodel_opt[modindex] << ".tif";
+ input=outputstream.str();
+ }
+ ImgReaderGdal imgReaderEst(input);
+
+ vector<double> obsBuffer;
+ vector<double> estReadBuffer;
+ vector<double> uncertReadBuffer;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
+ assert(irow<imgReaderEst.nrOfRow());
+ imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
+ imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
+ if(update){
+ imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow);
}
- else{//time update
+ for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
+ double kalmanGain=0;
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 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;
- 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);
+ 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];
}
- uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
}
+ imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
+ imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
+ progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
+ pfnProgress(progress,pszMessage,pProgressArg);
}
- imgWriterEst.writeData(estWriteBuffer,GDT_Float64,irow,0);
- imgWriterEst.writeData(uncertWriteBuffer,GDT_Float64,irow,1);
- progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
- pfnProgress(progress,pszMessage,pProgressArg);
- }
- imgWriterEst.close();
- imgReaderEst.close();
+ imgWriterEst.close();
+ imgReaderEst.close();
- if(update){
- imgReaderObs.close();
- ++obsindex;
+ if(update){
+ imgReaderObs.close();
+ --obsindex;
+ }
+ imgReaderModel1.close();
+ imgReaderModel2.close();
}
- 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