[pktools] 288/375: added support for uncertainty as second band in observation image for 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 179fdf25d103180f4be83fd9c657f0716e6e49be
Author: user <user at osgeolive.(none)>
Date: Mon Jun 16 21:45:06 2014 +0200
added support for uncertainty as second band in observation image for pkkalman.cc
---
doc/logo.odg | Bin 13276 -> 12241 bytes
src/algorithms/Filter.cc | 36 +++++++++++++++++
src/algorithms/Filter.h | 4 +-
src/algorithms/Filter2d.h | 3 +-
src/apps/pkfilter.cc | 19 ++++++++-
src/apps/pkkalman.cc | 98 ++++++++++++++++++++++++++++++++++------------
6 files changed, 130 insertions(+), 30 deletions(-)
diff --git a/doc/logo.odg b/doc/logo.odg
index 16b0427..55d33d4 100644
Binary files a/doc/logo.odg and b/doc/logo.odg differ
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 6644631..4095566 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -142,6 +142,42 @@ void filter::Filter::dwtCut(const ImgReaderGdal& input, ImgWriterGdal& output, c
}
}
+void filter::Filter::dwtCutFrom(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, int band){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+ vector<double> pixelInput(input.nrOfBand());
+ for(int x=0;x<input.nrOfCol();++x){
+ pixelInput=lineInput.selectCol(x);
+ dwtForward(pixelInput,wavelet_type,family);
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ if(iband>=band)
+ pixelInput[iband]=0;
+ }
+ dwtInverse(pixelInput,wavelet_type,family);
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ lineOutput[iband][x]=pixelInput[iband];
+ }
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ try{
+ output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in band " << iband << ", line " << y << endl;
+ }
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+}
+
void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
int origsize=data.size();
//make sure data size if power of 2
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 3cdab8e..38ce28e 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -33,7 +33,7 @@ extern "C" {
namespace filter
{
- enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28};
+ enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28, dwt_cut_from=29};
class Filter
{
@@ -77,6 +77,7 @@ public:
void dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family);
void dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family);
void dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut);
+ void dwtCutFrom(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, int band);
private:
@@ -85,6 +86,7 @@ private:
m_filterMap["dwt"]=filter::dwt;
m_filterMap["dwti"]=filter::dwti;
m_filterMap["dwt_cut"]=filter::dwt_cut;
+ m_filterMap["dwt_cut_from"]=filter::dwt_cut_from;
m_filterMap["stdev"]=filter::stdev;
m_filterMap["var"]=filter::var;
m_filterMap["min"]=filter::min;
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index dec5a32..deb27d8 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -58,7 +58,7 @@ extern "C" {
namespace filter2d
{
- enum FILTER_TYPE { median=100, var=101 , min=102, max=103, sum=104, mean=105, minmax=106, dilate=107, erode=108, close=109, open=110, homog=111, sobelx=112, sobely=113, sobelxy=114, sobelyx=115, smooth=116, density=117, majority=118, mixed=119, threshold=120, ismin=121, ismax=122, heterog=123, order=124, stdev=125, mrf=126, dwt=127, dwti=128, dwt_cut=129, scramble=130, shift=131, linearfeature=132, smoothnodata=133, countid=134};
+ enum FILTER_TYPE { median=100, var=101 , min=102, max=103, sum=104, mean=105, minmax=106, dilate=107, erode=108, close=109, open=110, homog=111, sobelx=112, sobely=113, sobelxy=114, sobelyx=115, smooth=116, density=117, majority=118, mixed=119, threshold=120, ismin=121, ismax=122, heterog=123, order=124, stdev=125, mrf=126, dwt=127, dwti=128, dwt_cut=129, scramble=130, shift=131, linearfeature=132, smoothnodata=133, countid=134, dwt_cut_from=135};
enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };//bicubic not supported yet...
@@ -162,6 +162,7 @@ private:
m_filterMap["dwt"]=filter2d::dwt;
m_filterMap["dwti"]=filter2d::dwti;
m_filterMap["dwt_cut"]=filter2d::dwt_cut;
+ m_filterMap["dwt_cut_from"]=filter2d::dwt_cut_from;
m_filterMap["scramble"]=filter2d::scramble;
m_filterMap["shift"]=filter2d::shift;
m_filterMap["linearfeature"]=filter2d::linearfeature;
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 8d0a63f..4ea7e06 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -42,7 +42,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
// Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation (North=0, East=90, South=180, West=270).");
- Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, countid, majority voting (only for classes), smoothnodata (smooth nodata values only) values, threshold local [...]
+ Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, countid, majority voting (only for classes), smoothnodata (smooth nodata values only) values, threshold local [...]
Optionpk<std::string> resample_opt("r", "resampling-method", "Resampling method for shifting operation (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<double> dimX_opt("dx", "dx", "filter kernel size in x, better use odd value to avoid image shift", 3);
Optionpk<double> dimY_opt("dy", "dy", "filter kernel size in y, better use odd value to avoid image shift", 3);
@@ -50,7 +50,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment, see also http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
Optionpk<short> class_opt("class", "class", "class value(s) to use for density, erosion, dilation, openening and closing, thresholding");
- Optionpk<double> threshold_opt("t", "threshold", "threshold value(s) to use for threshold filter (one for each class), or threshold to cut for dwt_cut (use 0 to keep all), or sigma for shift", 0);
+ Optionpk<double> threshold_opt("t", "threshold", "threshold value(s) to use for threshold filter (one for each class), or threshold to cut for dwt_cut (use 0 to keep all) or dwt_cut_from, or sigma for shift", 0);
Optionpk<short> nodata_opt("nodata", "nodata", "nodata value(s) for smoothnodata filter");
Optionpk<std::string> tap_opt("tap", "tap", "text file containing taps used for spatial filtering (from ul to lr). Use dimX and dimY to specify tap dimensions in x and y. Leave empty for not using taps");
Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
@@ -684,6 +684,21 @@ int main(int argc,char **argv) {
else
filter2d.dwtCut(input, output, wavelet_type_opt[0], family_opt[0], threshold_opt[0]);
break;
+ case(filter2d::dwt_cut_from):
+ if(down_opt[0]!=1){
+ std::cerr << "Error: down option not supported for this filter" << std::endl;
+ exit(1);
+ }
+ if(dimZ_opt.size()){
+ if(verbose_opt[0])
+ std::cout<< "DWT approximation in spectral domain" << std::endl;
+ filter1d.dwtCutFrom(input, output, wavelet_type_opt[0], family_opt[0], static_cast<int>(threshold_opt[0]));
+ }
+ else{
+ std::cerr << "Error: this filter is not supported in 2D" << std::endl;
+ exit(1);
+ }
+ break;
case(filter2d::threshold):
filter2d.setThresholds(threshold_opt);//deliberate fall through
case(filter2d::density):
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index daa6a5b..dd666c2 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -28,6 +28,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "algorithms/ImgRegression.h"
//todo: keep original resolution of coarse model raster dataset
+ //interprete 2nd band of obs dataset as uncert
using namespace std;
/*------------------
Main procedure
@@ -38,6 +39,7 @@ int main(int argc,char **argv) {
Optionpk<string> observation_opt("obs","observation","observation input datasets, e.g., landsat (use: -obs obs1 -obs obs2 etc.");
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> projection_opt("a_srs", "a_srs", "Override the projection for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid");
Optionpk<string> outputfw_opt("ofw", "outputfw", "Output raster dataset for forward model");
Optionpk<string> outputbw_opt("obw", "outputbw", "Output raster dataset for backward model");
Optionpk<string> outputfb_opt("ofb", "outputfb", "Output raster dataset for smooth model");
@@ -62,6 +64,7 @@ int main(int argc,char **argv) {
observation_opt.retrieveOption(argc,argv);
tmodel_opt.retrieveOption(argc,argv);
tobservation_opt.retrieveOption(argc,argv);
+ projection_opt.retrieveOption(argc,argv);
outputfw_opt.retrieveOption(argc,argv);
outputbw_opt.retrieveOption(argc,argv);
outputfb_opt.retrieveOption(argc,argv);
@@ -153,8 +156,13 @@ int main(int argc,char **argv) {
ImgWriterGdal imgWriterEst;
imgReaderObs.open(observation_opt[0]);
+
int ncol=imgReaderObs.nrOfCol();
int nrow=imgReaderObs.nrOfRow();
+ if(projection_opt.empty())
+ projection_opt.push_back(imgReaderObs.getProjection());
+ double geotransform[6];
+ imgReaderObs.getGeoTransform(geotransform);
string imageType=imgReaderObs.getImageType();
if(oformat_opt.size())//default
@@ -212,6 +220,8 @@ int main(int argc,char **argv) {
if(verbose_opt[0])
cout << "Opening image " << output << " for writing " << endl;
imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+ imgWriterEst.setProjectionProj4(projection_opt[0]);
+ imgWriterEst.setGeoTransform(geotransform);
imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
if(verbose_opt[0]){
@@ -243,27 +253,34 @@ int main(int argc,char **argv) {
}
else{//we have an observation at time 0
imgReaderObs.open(observation_opt[0]);
+ imgReaderObs.getGeoTransform(geotransform);
imgReaderObs.setNoData(obsnodata_opt);
for(int irow=0;irow<nrow;++irow){
vector<double> estReadBuffer;
imgReaderModel1.readData(estReadBuffer,GDT_Float64,irow);
vector<double> estWriteBuffer(ncol);
vector<double> uncertWriteBuffer(ncol);
+ vector<double> uncertObsBuffer;
imgReaderObs.readData(estWriteBuffer,GDT_Float64,irow,0);
+ if(imgReaderObs.nrOfBand()>1)
+ imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
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;
+ 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_opt[0]*uncertObs_opt[0]/noemer;
+ estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;
}
else{
//no need to fill write buffer (already done in imgReaderObs.readData
- uncertWriteBuffer[icol]=uncertObs_opt[0];
+ uncertWriteBuffer[icol]=uncertObs;
}
}
}
@@ -296,6 +313,9 @@ int main(int argc,char **argv) {
//two band output band0=estimation, band1=uncertainty
imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+ imgWriterEst.setProjectionProj4(projection_opt[0]);
+ imgWriterEst.setGeoTransform(geotransform);
+
imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
//calculate regression between two subsequence model inputs
@@ -320,6 +340,7 @@ int main(int argc,char **argv) {
cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.getGeoTransform(geotransform);
imgReaderObs.setNoData(obsnodata_opt);
//calculate regression between model and observation
errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
@@ -335,6 +356,7 @@ int main(int argc,char **argv) {
ImgReaderGdal imgReaderEst(input);
vector<double> obsBuffer;
+ vector<double> uncertObsBuffer;
vector<double> estReadBuffer;
vector<double> uncertReadBuffer;
vector<double> estWriteBuffer(ncol);
@@ -344,7 +366,9 @@ int main(int argc,char **argv) {
imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
if(update){
- imgReaderObs.readData(obsBuffer,GDT_Float64,irow);
+ imgReaderObs.readData(obsBuffer,GDT_Float64,irow,0);
+ if(imgReaderObs.nrOfBand()>1)
+ imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
}
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
double estValue=estReadBuffer[icol];
@@ -370,8 +394,11 @@ int main(int argc,char **argv) {
//observation update
if(update&&!imgReaderObs.isNoData(obsBuffer[icol])){
double kalmanGain=1;
- if((uncertWriteBuffer[icol]+uncertObs_opt[0])>eps_opt[0])
- kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs_opt[0]);
+ double uncertObs=uncertObs_opt[0];
+ if(uncertObsBuffer.size()>icol)
+ uncertObs=uncertObsBuffer[icol];
+ if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0])
+ kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs);
assert(kalmanGain<=1);
estWriteBuffer[icol]+=kalmanGain*(obsBuffer[icol]-estWriteBuffer[icol]);
uncertWriteBuffer[icol]*=(1-kalmanGain);
@@ -411,6 +438,8 @@ int main(int argc,char **argv) {
if(verbose_opt[0])
cout << "Opening image " << output << " for writing " << endl;
imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+ imgWriterEst.setProjectionProj4(projection_opt[0]);
+ imgWriterEst.setGeoTransform(geotransform);
imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
if(verbose_opt[0]){
@@ -441,27 +470,34 @@ int main(int argc,char **argv) {
}
else{//we have an observation at end time
imgReaderObs.open(observation_opt.back());
+ imgReaderObs.getGeoTransform(geotransform);
imgReaderObs.setNoData(obsnodata_opt);
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(estReadBuffer,GDT_Float64,irow);
+ vector<double> uncertObsBuffer;
+ imgReaderObs.readData(estReadBuffer,GDT_Float64,irow,0);
+ if(imgReaderObs.nrOfBand()>1)
+ imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
for(int icol=0;icol<imgWriterEst.nrOfCol();++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;
+ 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_opt[0]*uncertObs_opt[0]/noemer;
+ estWriteBuffer[icol]+=uncertModel_opt[0]*uncertObs*uncertObs/noemer;
}
else{
//no need to fill write buffer (already done in imgReaderObs.readData
- uncertWriteBuffer[icol]=uncertObs_opt[0];
+ uncertWriteBuffer[icol]=uncertObs;
}
}
}
@@ -493,6 +529,8 @@ int main(int argc,char **argv) {
}
//two band output band0=estimation, band1=uncertainty
imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+ imgWriterEst.setProjectionProj4(projection_opt[0]);
+ imgWriterEst.setGeoTransform(geotransform);
imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
//calculate regression between two subsequence model inputs
@@ -516,6 +554,7 @@ int main(int argc,char **argv) {
if(verbose_opt[0])
cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.getGeoTransform(geotransform);
imgReaderObs.setNoData(obsnodata_opt);
//calculate regression between model and observation
errObs=imgreg.getRMSE(imgReaderModel1,imgReaderObs,c0obs,c1obs,verbose_opt[0]);
@@ -533,6 +572,7 @@ int main(int argc,char **argv) {
ImgReaderGdal imgReaderEst(input);
vector<double> obsBuffer;
+ vector<double> uncertObsBuffer;
vector<double> estReadBuffer;
vector<double> uncertReadBuffer;
vector<double> estWriteBuffer(ncol);
@@ -543,7 +583,9 @@ int main(int argc,char **argv) {
imgReaderEst.readData(estReadBuffer,GDT_Float64,irow,0);
imgReaderEst.readData(uncertReadBuffer,GDT_Float64,irow,1);
if(update){
- imgReaderObs.readData(obsBuffer,GDT_Float64,irow);
+ imgReaderObs.readData(obsBuffer,GDT_Float64,irow,0);
+ if(imgReaderObs.nrOfBand()>1)
+ imgReaderObs.readData(uncertObsBuffer,GDT_Float64,irow,1);
}
for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
double estValue=estReadBuffer[icol];
@@ -560,17 +602,20 @@ int main(int argc,char **argv) {
if(errMod<eps_opt[0])
totalUncertainty=errObs;
else if(errObs<eps_opt[0])
- totalUncertainty=errObs;
+ totalUncertainty=errObs;
else{
- totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
- totalUncertainty=sqrt(1.0/totalUncertainty);
- }
+ totalUncertainty=1.0/errMod/errMod+1/errObs/errObs;
+ totalUncertainty=sqrt(1.0/totalUncertainty);
+ }
uncertWriteBuffer[icol]=totalUncertainty+uncertReadBuffer[icol];
//observation update
if(update&&!imgReaderObs.isNoData(obsBuffer[icol])){
double kalmanGain=1;
- if((uncertWriteBuffer[icol]+uncertObs_opt[0])>eps_opt[0])
- kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs_opt[0]);
+ double uncertObs=uncertObs_opt[0];
+ if(uncertObsBuffer.size()>icol)
+ uncertObs=uncertObsBuffer[icol];
+ if((uncertWriteBuffer[icol]+uncertObs)>eps_opt[0])
+ kalmanGain=uncertWriteBuffer[icol]/(uncertWriteBuffer[icol]+uncertObs);
assert(kalmanGain<=1);
estWriteBuffer[icol]+=kalmanGain*(obsBuffer[icol]-estWriteBuffer[icol]);
uncertWriteBuffer[icol]*=(1-kalmanGain);
@@ -617,6 +662,8 @@ int main(int argc,char **argv) {
//two band output band0=estimation, band1=uncertainty
imgWriterEst.open(output,ncol,nrow,2,GDT_Float32,imageType,option_opt);
+ imgWriterEst.setProjectionProj4(projection_opt[0]);
+ imgWriterEst.setGeoTransform(geotransform);
imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
//open forward and backward estimates
@@ -660,6 +707,7 @@ int main(int argc,char **argv) {
if(verbose_opt[0])
cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.getGeoTransform(geotransform);
imgReaderObs.setNoData(obsnodata_opt);
//calculate regression between model and observation
}
@@ -685,22 +733,20 @@ int main(int argc,char **argv) {
double B=estBackwardBuffer[icol];
double C=uncertForwardBuffer[icol]*uncertForwardBuffer[icol];
double D=uncertBackwardBuffer[icol]*uncertBackwardBuffer[icol];
- double uncertObs=0;
+ double uncertObs=uncertObs_opt[0];
if(update){//check for nodata in observation
- if(imgReaderObs.nrOfBand()>1)
- uncertObs=uncertObsBuffer[icol];
- else if(imgReaderObs.isNoData(estWriteBuffer[icol]))
+ if(imgReaderObs.isNoData(estWriteBuffer[icol]))
uncertObs=uncertNodata_opt[0];
- else
- uncertObs=uncertObs_opt[0];
+ else if(uncertObsBuffer.size()>icol)
+ uncertObs=uncertObsBuffer[icol];
}
double noemer=(C+D);
//todo: consistently check for division by zero...
if(noemer<eps_opt[0]){//simple average if both uncertainties are ~>0
estWriteBuffer[icol]=0.5*(A+B);
- uncertWriteBuffer[icol]=uncertObs_opt[0];
+ uncertWriteBuffer[icol]=uncertObs;
}
else{
estWriteBuffer[icol]=(A*D+B*C)/noemer;
--
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