pktools  2.6.7
Processing Kernel for geospatial data
pkkalman.cc
1 /**********************************************************************
2 pkkalman.cc: produce kalman filtered raster time series
3 Copyright (C) 2008-2014 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include <sstream>
21 #include <vector>
22 #include <algorithm>
23 #include "base/Optionpk.h"
24 #include "base/Vector2d.h"
25 #include "imageclasses/ImgReaderGdal.h"
26 #include "imageclasses/ImgWriterGdal.h"
27 #include "imageclasses/ImgUpdaterGdal.h"
28 #include "algorithms/StatFactory.h"
29 
30 /******************************************************************************/
91 using namespace std;
92 /*------------------
93  Main procedure
94  ----------------*/
95 int main(int argc,char **argv) {
96  Optionpk<string> direction_opt("dir","direction","direction to run model (forward|backward|smooth)","forward");
97  Optionpk<string> model_opt("mod","model","coarse spatial resolution input datasets(s) used as model. Use either multi-band input (-model multiband_model.tif) or multiple single-band inputs (-mod model1 -mod model2 etc.)");
98  Optionpk<string> modelmask_opt("modmask","modmask","model mask datasets(s). Must have same dimension as model input. Use either multi-band input or multiple single-band inputs");
99  Optionpk<string> observation_opt("obs","observation","fine spatial resolution input dataset(s) used as observation. Use either multi-band input (-obs multiband_obs.tif) or multiple single-band inputs (-obs obs1 -obs obs2 etc.)");
100  Optionpk<string> observationmask_opt("obsmask","obsmask","observation mask dataset(s). Must have same dimension as observation input (use multi-band input or multiple single-band inputs");
101  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.");
102  Optionpk<int> tobservation_opt("tobs","tobservation","time sequence of observation input. Sequence must have exact same length as observation input)");
103  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");
104  Optionpk<string> outputfw_opt("ofw", "outputfw", "Output raster dataset for forward model");
105  Optionpk<string> uncertfw_opt("u_ofw", "u_outputfw", "Uncertainty output raster dataset for forward model");
106  Optionpk<string> outputbw_opt("obw", "outputbw", "Output raster dataset for backward model");
107  Optionpk<string> uncertbw_opt("u_obw", "u_outputbw", "Uncertainty output raster dataset for backward model");
108  Optionpk<string> outputfb_opt("ofb", "outputfb", "Output raster dataset for smooth model");
109  Optionpk<string> uncertfb_opt("u_ofb", "u_outputfb", "Uncertainty output raster dataset for smooth model");
110  Optionpk<string> gain_opt("gain", "gain", "Output raster dataset for gain");
111  Optionpk<double> modnodata_opt("modnodata", "modnodata", "invalid value for model input", 0);
112  Optionpk<double> obsnodata_opt("obsnodata", "obsnodata", "invalid value for observation input", 0);
113  Optionpk<double> obsmin_opt("obsmin", "obsmin", "Minimum value for observation data");
114  Optionpk<double> obsmax_opt("obsmax", "obsmax", "Maximum value for observation data");
115  Optionpk<double> msknodata_opt("msknodata", "msknodata", "Mask value not to consider", 0);
116  Optionpk<short> mskband_opt("mskband", "mskband", "Mask band to read (0 indexed)", 0);
117  Optionpk<double> eps_opt("eps", "eps", "epsilon for non zero division", 0.00001);
118  Optionpk<double> uncertModel_opt("um", "uncertmodel", "Uncertainty of model",1);
119  Optionpk<double> uncertObs_opt("uo", "uncertobs", "Uncertainty of valid observations",1);
120  Optionpk<double> processNoise_opt("q", "q", "Process noise: expresses instability (variance) of proportions of fine res pixels within a moderate resolution pixel",1);
121  Optionpk<double> uncertNodata_opt("unodata", "uncertnodata", "Uncertainty in case of no-data values in observation", 100);
122  Optionpk<int> down_opt("down", "down", "Downsampling factor for reading model data to calculate regression");
123  Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
124  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff",2);
125  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
126  Optionpk<short> verbose_opt("v", "verbose", "verbose mode when positive", 0);
127 
128  observationmask_opt.setHide(1);
129  modelmask_opt.setHide(1);
130  tmodel_opt.setHide(1);
131  tobservation_opt.setHide(1);
132  projection_opt.setHide(1);
133  uncertfw_opt.setHide(1);
134  uncertbw_opt.setHide(1);
135  uncertfb_opt.setHide(1);
136  obsmin_opt.setHide(1);
137  obsmax_opt.setHide(1);
138  msknodata_opt.setHide(1);
139  mskband_opt.setHide(1);
140  eps_opt.setHide(1);
141  uncertNodata_opt.setHide(1);
142  down_opt.setHide(1);
143  otype_opt.setHide(1);
144  oformat_opt.setHide(1);
145  option_opt.setHide(1);
146  verbose_opt.setHide(1);
147  gain_opt.setHide(2);
148 
149  bool doProcess;//stop process when program was invoked with help option (-h --help)
150  try{
151  doProcess=direction_opt.retrieveOption(argc,argv);
152  model_opt.retrieveOption(argc,argv);
153  modelmask_opt.retrieveOption(argc,argv);
154  observation_opt.retrieveOption(argc,argv);
155  observationmask_opt.retrieveOption(argc,argv);
156  tmodel_opt.retrieveOption(argc,argv);
157  tobservation_opt.retrieveOption(argc,argv);
158  projection_opt.retrieveOption(argc,argv);
159  outputfw_opt.retrieveOption(argc,argv);
160  uncertfw_opt.retrieveOption(argc,argv);
161  outputbw_opt.retrieveOption(argc,argv);
162  uncertbw_opt.retrieveOption(argc,argv);
163  outputfb_opt.retrieveOption(argc,argv);
164  uncertfb_opt.retrieveOption(argc,argv);
165  gain_opt.retrieveOption(argc,argv);
166  modnodata_opt.retrieveOption(argc,argv);
167  obsnodata_opt.retrieveOption(argc,argv);
168  obsmin_opt.retrieveOption(argc,argv);
169  obsmax_opt.retrieveOption(argc,argv);
170  msknodata_opt.retrieveOption(argc,argv);
171  mskband_opt.retrieveOption(argc,argv);
172  eps_opt.retrieveOption(argc,argv);
173  uncertModel_opt.retrieveOption(argc,argv);
174  uncertObs_opt.retrieveOption(argc,argv);
175  processNoise_opt.retrieveOption(argc,argv);
176  uncertNodata_opt.retrieveOption(argc,argv);
177  down_opt.retrieveOption(argc,argv);
178  otype_opt.retrieveOption(argc,argv);
179  oformat_opt.retrieveOption(argc,argv);
180  option_opt.retrieveOption(argc,argv);
181  verbose_opt.retrieveOption(argc,argv);
182  }
183  catch(string predefinedString){
184  std::cout << predefinedString << std::endl;
185  exit(0);
186  }
187  if(!doProcess){
188  std::cerr << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
189  exit(0);//help was invoked, stop processing
190  }
191 
192  try{
193  ostringstream errorStream;
194  if(model_opt.size()<1){
195  errorStream << "Error: no model dataset selected, use option -mod" << endl;
196  throw(errorStream.str());
197  }
198  if(observation_opt.size()<1){
199  errorStream << "Error: no observation dataset selected, use option -obs" << endl;
200  throw(errorStream.str());
201  }
202  if(find(direction_opt.begin(),direction_opt.end(),"forward")!=direction_opt.end()){
203  if(outputfw_opt.empty()){
204  errorStream << "Error: output forward datasets is not provided, use option -ofw" << endl;
205  throw(errorStream.str());
206  }
207  if(uncertfw_opt.empty()){
208  ostringstream uncertStream;
209  uncertStream << outputfw_opt[0] << "_uncert";
210  uncertfw_opt.push_back(uncertStream.str());
211  }
212  }
213  if(find(direction_opt.begin(),direction_opt.end(),"backward")!=direction_opt.end()){
214  if(outputbw_opt.empty()){
215  errorStream << "Error: output backward datasets is not provided, use option -obw" << endl;
216  throw(errorStream.str());
217  }
218  if(uncertbw_opt.empty()){
219  ostringstream uncertStream;
220  uncertStream << outputbw_opt[0] << "_uncert";
221  uncertbw_opt.push_back(uncertStream.str());
222  }
223  }
224  // if(model_opt.size()<observation_opt.size()){
225  // errorStream << "Error: sequence of models should be larger than observations" << endl;
226  // throw(errorStream.str());
227  // }
228  if(tmodel_opt.empty()){
229  cout << "Warning: model time sequence is not provided, self generating time sequence from model input" << endl;
230  }
231  if(tobservation_opt.empty()){
232  cout << "Warning: observation time sequence is not provided, self generating time sequence from observation input" << endl;
233  }
234  if(find(direction_opt.begin(),direction_opt.end(),"smooth")!=direction_opt.end()){
235  if(outputfw_opt.empty()){
236  errorStream << "Error: output forward dataset is not provided, use option -ofw" << endl;
237  throw(errorStream.str());
238  }
239  if(outputbw_opt.empty()){
240  errorStream << "Error: output backward datasets is not provided, use option -obw" << endl;
241  throw(errorStream.str());
242  }
243  if(outputfb_opt.empty()){
244  errorStream << "Error: output smooth datasets is not provided, use option -ofb" << endl;
245  throw(errorStream.str());
246  }
247  if(uncertfb_opt.empty()){
248  ostringstream uncertStream;
249  uncertStream << outputfb_opt[0] << "_uncert";
250  uncertfb_opt.push_back(uncertStream.str());
251  }
252  }
253  }
254  catch(string errorString){
255  std::cout << errorString << std::endl;
256  exit(1);
257  }
258 
260  stat.setNoDataValues(modnodata_opt);
261  ImgReaderGdal imgReaderModel1;
262  ImgReaderGdal imgReaderModel2;
263  ImgReaderGdal imgReaderModel1Mask;
264  ImgReaderGdal imgReaderModel2Mask;
265  ImgReaderGdal imgReaderObs;
266  ImgReaderGdal imgReaderObsMask;
267  //test
268  ImgWriterGdal imgWriterGain;
269 
270  imgReaderModel1.open(model_opt[0]);
271  imgReaderModel1.setNoData(modnodata_opt);
272  imgReaderObs.open(observation_opt[0]);
273  imgReaderObs.setNoData(obsnodata_opt);
274  // if(observationmask_opt.empty())
275  // observationmask_opt=observation_opt;
276  if(modelmask_opt.size()){
277  imgReaderModel1Mask.open(modelmask_opt[0]);
278  imgReaderModel1Mask.setNoData(msknodata_opt);
279  }
280  if(observationmask_opt.size()){
281  imgReaderObsMask.open(observationmask_opt[0]);
282  imgReaderObsMask.setNoData(msknodata_opt);
283  }
284 
285  unsigned int nobs=(observation_opt.size()>1)? observation_opt.size() : imgReaderObs.nrOfBand();
286  unsigned int nmodel=(model_opt.size()>1)? model_opt.size() : imgReaderModel1.nrOfBand();
287 
288  if(verbose_opt[0]){
289  cout << "number of observations: " << nobs << endl;
290  cout << "number of models: " << nmodel << endl;
291  }
292 
293  int ncol=imgReaderObs.nrOfCol();
294  int nrow=imgReaderObs.nrOfRow();
295  if(projection_opt.empty())
296  projection_opt.push_back(imgReaderObs.getProjection());
297  double geotransform[6];
298  imgReaderObs.getGeoTransform(geotransform);
299 
300  GDALDataType theType=GDT_Unknown;
301  if(verbose_opt[0])
302  cout << "possible output data types: ";
303  for(int iType = 0; iType < GDT_TypeCount; ++iType){
304  if(verbose_opt[0])
305  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
306  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
307  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
308  otype_opt[0].c_str()))
309  theType=(GDALDataType) iType;
310  }
311  if(theType==GDT_Unknown)
312  theType=imgReaderObs.getDataType();
313 
314  string imageType;//=imgReaderObs.getImageType();
315  if(oformat_opt.size())//default
316  imageType=oformat_opt[0];
317  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
318  string theInterleave="INTERLEAVE=";
319  theInterleave+=imgReaderObs.getInterleave();
320  option_opt.push_back(theInterleave);
321  }
322 
323  if(down_opt.empty()){
324  double resModel=imgReaderModel1.getDeltaX();
325  double resObs=imgReaderObs.getDeltaX();
326  int down=static_cast<int>(ceil(resModel/resObs));
327  if(!(down%2))
328  down+=1;
329  down_opt.push_back(down);
330  }
331 
332  int obsindex=0;
333 
334  const char* pszMessage;
335  void* pProgressArg=NULL;
336  GDALProgressFunc pfnProgress=GDALTermProgress;
337  double progress=0;
338 
339  double errObs=uncertNodata_opt[0];//start with high initial value in case we do not have first observation at time 0
340 
341  while(tmodel_opt.size()<nmodel)
342  tmodel_opt.push_back(tmodel_opt.size());
343  try{
344  if(tobservation_opt.size()<nobs){
345  if(nobs==nmodel){
346  while(tobservation_opt.size()<nobs)
347  tobservation_opt.push_back(tobservation_opt.size());
348  }
349  else{
350  ostringstream errorStream;
351  errorStream << "Error: please provide time sequence for observation using option tobs" << endl;
352  throw(errorStream.str());
353  }
354  }
355  }
356  catch(string errorString){
357  std::cout << errorString << std::endl;
358  exit(1);
359  }
360 
361  vector<int> relobsindex;
362 
363  for(int tindex=0;tindex<tobservation_opt.size();++tindex){
364  vector<int>::iterator modit;
365  modit=upper_bound(tmodel_opt.begin(),tmodel_opt.end(),tobservation_opt[tindex]);
366  int relpos=modit-tmodel_opt.begin()-1;
367  assert(relpos>=0);//todo: for now, we assume model is available at time before first measurement
368  relobsindex.push_back(relpos);
369  if(verbose_opt[0]){
370  cout << "observation " << tindex << ": " << "relative position in model time series is " << relpos << ", date of observation is (tobservation_opt[tindex]): " << tobservation_opt[tindex] << ", relobsindex.back(): " << relobsindex.back();
371  if(observation_opt.size()>tindex)
372  cout << ", filename observation: " << observation_opt[tindex];
373  else
374  cout << ", observation band index: " << tindex;
375  if(model_opt.size()>relpos)
376  cout << ", filename of corresponding model: " << model_opt[relpos] << endl;
377  else
378  cout << ", band index of corresponding model: " << relpos;
379  }
380  }
381 
382  int ndigit=log(1.0*tmodel_opt.back())/log(10.0)+1;
383 
384  double geox=0;
385  double geoy=0;
386 
387  if(model_opt.size()==nmodel)
388  imgReaderModel1.close();
389  if(modelmask_opt.size()==nmodel)
390  imgReaderModel1Mask.close();
391  if(observation_opt.size()==nobs)
392  imgReaderObs.close();
393  if(observationmask_opt.size()==nobs)
394  imgReaderObsMask.close();
395 
396  try{
397  if(find(direction_opt.begin(),direction_opt.end(),"forward")!=direction_opt.end()){
399  cout << "Running forward model" << endl;
400  obsindex=0;
401  if(verbose_opt[0])
402  cout << "Opening image " << outputfw_opt[0] << " for writing " << endl << flush;
403 
404  ImgWriterGdal imgWriterEst;
405  ImgWriterGdal imgWriterUncert;
406  imgWriterEst.open(outputfw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
407  imgWriterEst.setProjectionProj4(projection_opt[0]);
408  imgWriterEst.setGeoTransform(geotransform);
409  imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
410  imgWriterUncert.open(uncertfw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
411  imgWriterUncert.setProjectionProj4(projection_opt[0]);
412  imgWriterUncert.setGeoTransform(geotransform);
413 
414  try{
415  //test
416  if(gain_opt.size()){
417  imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
418  imgWriterGain.setProjectionProj4(projection_opt[0]);
419  imgWriterGain.setGeoTransform(geotransform);
420  imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
421  }
422 
423  if(verbose_opt[0]){
424  cout << "processing time " << tmodel_opt[0] << endl;
425  if(obsindex<relobsindex.size()){
426  assert(tmodel_opt.size()>relobsindex[obsindex]);
427  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
428  }
429  else
430  cout << "There is no next observation" << endl;
431  }
432  if(model_opt.size()==nmodel){
433  imgReaderModel1.open(model_opt[0]);
434  imgReaderModel1.setNoData(modnodata_opt);
435  }
436  if(modelmask_opt.size()==nmodel){
437  imgReaderModel1Mask.open(modelmask_opt[0]);
438  imgReaderModel1Mask.setNoData(msknodata_opt);
439  }
440  }
441  catch(string errorString){
442  cerr << errorString << endl;
443  }
444  catch(...){
445  cerr << "Error opening file " << model_opt[0] << endl;
446  }
447 
448  double modRow=0;
449  double modCol=0;
450  double lowerCol=0;
451  double upperCol=0;
452  RESAMPLE theResample=BILINEAR;
453 
454  if(relobsindex[0]>0){//initialize output_opt[0] as model[0]
455  //write first model as output
456  if(verbose_opt[0])
457  cout << "write first model as output" << endl;
458  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
459  vector<double> estReadBuffer;
460  vector<double> lineModelMask;
461  vector<double> estWriteBuffer(ncol);
462  vector<double> uncertWriteBuffer(ncol);
463  //test
464  vector<double> gainWriteBuffer(ncol);
465  try{
466  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
467  imgWriterEst.image2geo(0,irow,geox,geoy);
468  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
469  if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
470  cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
471  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
472  }
473 
474  int readModelBand=(model_opt.size()==nmodel)? 0:0;
475  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
476  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
477  if(modelmask_opt.size())
478  imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
479  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
480  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
481  imgWriterEst.image2geo(icol,irow,geox,geoy);
482  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
483  if(modelmask_opt.size()){
484  if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
485  estWriteBuffer[icol]=obsnodata_opt[0];
486  uncertWriteBuffer[icol]=uncertNodata_opt[0];
487  //test
488  gainWriteBuffer[icol]=obsnodata_opt[0];
489  continue;
490  }
491  }
492  lowerCol=modCol-0.5;
493  lowerCol=static_cast<int>(lowerCol);
494  upperCol=modCol+0.5;
495  upperCol=static_cast<int>(upperCol);
496  if(lowerCol<0)
497  lowerCol=0;
498  if(upperCol>=imgReaderModel1.nrOfCol())
499  upperCol=imgReaderModel1.nrOfCol()-1;
500  double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
501  // double modValue=estReadBuffer[modCol];
502  if(imgReaderModel1.isNoData(modValue)){
503  estWriteBuffer[icol]=obsnodata_opt[0];
504  uncertWriteBuffer[icol]=uncertNodata_opt[0];
505  //test
506  gainWriteBuffer[icol]=obsnodata_opt[0];
507  continue;
508  }
509  estWriteBuffer[icol]=modValue;
510  if(obsmin_opt.size()){
511  if(estWriteBuffer[icol]<obsmin_opt[0])
512  estWriteBuffer[icol]=obsmin_opt[0];
513  }
514  if(obsmax_opt.size()){
515  if(estWriteBuffer[icol]>obsmax_opt[0])
516  estWriteBuffer[icol]=obsmax_opt[0];
517  }
518  uncertWriteBuffer[icol]=uncertModel_opt[0];
519  //test
520  gainWriteBuffer[icol]=0;
521  }
522  }
523  imgWriterEst.writeData(estWriteBuffer,irow,0);
524  imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
525  //test
526  if(gain_opt.size())
527  imgWriterGain.writeData(gainWriteBuffer,irow,0);
528  }
529  }
530  catch(string errorString){
531  cerr << errorString << endl;
532  }
533  catch(...){
534  cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
535  }
536  }
537  }
538  else{//we have a measurement
539  if(verbose_opt[0])
540  cout << "we have a measurement at initial time" << endl;
541  if(observation_opt.size()==nobs){
542  imgReaderObs.open(observation_opt[0]);
543  imgReaderObs.setNoData(obsnodata_opt);
544  }
545  if(observationmask_opt.size()==nobs){
546  imgReaderObsMask.open(observationmask_opt[0]);
547  imgReaderObsMask.setNoData(msknodata_opt);
548  }
549  imgReaderObs.getGeoTransform(geotransform);
550 
551  vector< vector<double> > obsLineVector(down_opt[0]);
552  vector<double> obsLineBuffer;
553  vector<double> obsMaskLineBuffer;
554  vector<double> modelMaskLineBuffer;
555  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
556  vector<double> estReadBuffer;
557  vector<double> estWriteBuffer(ncol);
558  vector<double> uncertWriteBuffer(ncol);
559  vector<double> uncertObsLineBuffer;
560  //test
561  vector<double> gainWriteBuffer(ncol);
562 
563  if(verbose_opt[0])
564  cout << "initialize obsLineVector" << endl;
565  assert(down_opt[0]%2);//window size must be odd
566  int readObsBand=(observation_opt.size()==nobs)? 0:0;
567  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:0;
568  int readModelBand=(model_opt.size()==nmodel)? 0:0;
569  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
570  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
571  if(iline<0)//replicate line 0
572  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
573  else
574  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
575  }
576  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
577  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
578  imgWriterEst.image2geo(0,irow,geox,geoy);
579  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
580  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
581  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
582  if(modelmask_opt.size())
583  imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
584  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
585  obsLineVector.erase(obsLineVector.begin());
586  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
587  obsLineVector.push_back(obsLineBuffer);
588 
589  if(observationmask_opt.size())
590  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
591 
592  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
593  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
594  imgWriterEst.image2geo(icol,irow,geox,geoy);
595  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
596  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
597  bool modelIsNoData=false;
598  if(modelmask_opt.size())
599  modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
600  lowerCol=modCol-0.5;
601  lowerCol=static_cast<int>(lowerCol);
602  upperCol=modCol+0.5;
603  upperCol=static_cast<int>(upperCol);
604  if(lowerCol<0)
605  lowerCol=0;
606  if(upperCol>=imgReaderModel1.nrOfCol())
607  upperCol=imgReaderModel1.nrOfCol()-1;
608  double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
609  // double modValue=estReadBuffer[modCol];
610  double errMod=uncertModel_opt[0];//*stdDev*stdDev;
611  modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
612  bool obsIsNoData=false;
613  if(observationmask_opt.size())
614  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
615  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
616  if(modelIsNoData){//model is nodata: retain observation
617  if(obsIsNoData){//both model and observation nodata
618  estWriteBuffer[icol]=obsnodata_opt[0];
619  uncertWriteBuffer[icol]=uncertNodata_opt[0];
620  //test
621  gainWriteBuffer[icol]=obsnodata_opt[0];
622  }
623  else{
624  estWriteBuffer[icol]=obsLineBuffer[icol];
625  if(obsmin_opt.size()){
626  if(estWriteBuffer[icol]<obsmin_opt[0])
627  estWriteBuffer[icol]=obsmin_opt[0];
628  }
629  if(obsmax_opt.size()){
630  if(estWriteBuffer[icol]>obsmax_opt[0])
631  estWriteBuffer[icol]=obsmax_opt[0];
632  }
633  uncertWriteBuffer[icol]=uncertObs_opt[0];
634  }
635  }
636  else{//model is valid: calculate estimate from model
637  estWriteBuffer[icol]=modValue;
638  uncertWriteBuffer[icol]=errMod;//in case observation is not valid
639  //test
640  gainWriteBuffer[icol]=0;
641  }
642  //measurement update
643  if(!obsIsNoData){
644  // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
645  double kalmanGain=1;
646  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
647  int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
648  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
649  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
650  obsWindowBuffer.clear();
651  for(int iline=0;iline<obsLineVector.size();++iline){
652  for(int isample=minCol;isample<=maxCol;++isample){
653  assert(isample<obsLineVector[iline].size());
654  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
655  }
656  }
657  if(!modelIsNoData){//model is valid
658  statfactory::StatFactory statobs;
659  statobs.setNoDataValues(obsnodata_opt);
660  double obsMeanValue=0;
661  double obsVarValue=0;
662  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
663  double difference=0;
664  difference=obsMeanValue-modValue;
665  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
666  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
667  // double errorCovariance=errMod;
668  double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
669  if(errorCovariance+errObs>eps_opt[0])
670  kalmanGain=errorCovariance/(errorCovariance+errObs);
671  else
672  kalmanGain=1;
673  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
674  if(obsmin_opt.size()){
675  if(estWriteBuffer[icol]<obsmin_opt[0])
676  estWriteBuffer[icol]=obsmin_opt[0];
677  }
678  if(obsmax_opt.size()){
679  if(estWriteBuffer[icol]>obsmax_opt[0])
680  estWriteBuffer[icol]=obsmax_opt[0];
681  if(uncertWriteBuffer[icol]>obsmax_opt[0])
682  uncertWriteBuffer[icol]=obsmax_opt[0];
683  }
684  }
685  if(kalmanGain>=1)
686  kalmanGain=1;
687  //test
688  gainWriteBuffer[icol]=kalmanGain;
689  }
690  }
691  }
692  imgWriterEst.writeData(estWriteBuffer,irow,0);
693  imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
694  //test
695  if(gain_opt.size())
696  imgWriterGain.writeData(gainWriteBuffer,irow,0);
697  }
698  }
699  if(observation_opt.size()==nobs)
700  imgReaderObs.close();
701  if(observationmask_opt.size()==nobs)
702  imgReaderObsMask.close();
703  ++obsindex;
704  }
705  if(model_opt.size()==nmodel)
706  imgReaderModel1.close();
707  if(modelmask_opt.size()==nmodel)
708  imgReaderModel1Mask.close();
709  imgWriterEst.close();
710  imgWriterUncert.close();
711 
712  ImgUpdaterGdal imgUpdaterEst;
713  ImgUpdaterGdal imgUpdaterUncert;
714  for(int modindex=1;modindex<nmodel;++modindex){
715  imgUpdaterEst.open(outputfw_opt[0]);
716  imgUpdaterEst.setNoData(obsnodata_opt);
717  imgUpdaterUncert.open(uncertfw_opt[0]);
718  if(verbose_opt[0]){
719  cout << "processing time " << tmodel_opt[modindex] << endl;
720  if(obsindex<relobsindex.size())
721  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
722  else
723  cout << "There is no next observation" << endl;
724  }
725 
726  //calculate regression between two subsequence model inputs
727  if(model_opt.size()==nmodel){
728  imgReaderModel1.open(model_opt[modindex-1]);
729  imgReaderModel1.setNoData(modnodata_opt);
730  imgReaderModel2.open(model_opt[modindex]);
731  imgReaderModel2.setNoData(modnodata_opt);
732  }
733  if(modelmask_opt.size()==nmodel){
734  imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
735  imgReaderModel1Mask.setNoData(msknodata_opt);
736  imgReaderModel2Mask.open(modelmask_opt[modindex]);
737  imgReaderModel2Mask.setNoData(msknodata_opt);
738  }
739 
740  pfnProgress(progress,pszMessage,pProgressArg);
741 
742  bool update=false;
743  if(obsindex<relobsindex.size()){
744  update=(relobsindex[obsindex]==modindex);
745  }
746  if(update){
747  if(observation_opt.size()==nobs){
748  if(verbose_opt[0])
749  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
750  imgReaderObs.open(observation_opt[obsindex]);
751  imgReaderObs.getGeoTransform(geotransform);
752  imgReaderObs.setNoData(obsnodata_opt);
753  }
754  if(observationmask_opt.size()==nobs){
755  imgReaderObsMask.open(observationmask_opt[obsindex]);
756  imgReaderObsMask.setNoData(msknodata_opt);
757  }
758  }
759  //prediction (also to fill cloudy pixels in measurement update mode)
760  string input;
761  input=outputfw_opt[0];
762 
763  vector< vector<double> > obsLineVector(down_opt[0]);
764  vector<double> obsLineBuffer;
765  vector<double> obsMaskLineBuffer;
766  vector<double> model1MaskLineBuffer;
767  vector<double> model2MaskLineBuffer;
768  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
769  vector<double> model1LineBuffer;
770  vector<double> model2LineBuffer;
771  vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
772  vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
773  vector<double> uncertObsLineBuffer;
774  vector< vector<double> > estLineVector(down_opt[0]);
775  vector<double> estLineBuffer;
776  vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
777  vector<double> uncertReadBuffer;
778  vector<double> estWriteBuffer(ncol);
779  vector<double> uncertWriteBuffer(ncol);
780  //test
781  vector<double> gainWriteBuffer(ncol);
782 
783  int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
784  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
785  int readModel1Band=(model_opt.size()==nmodel)? 0:modindex-1;
786  int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
787  int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex-1;
788  int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
789 
790  //initialize obsLineVector if update
791  if(update){
792  if(verbose_opt[0])
793  cout << "initialize obsLineVector" << endl;
794  assert(down_opt[0]%2);//window size must be odd
795  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
796  if(iline<0)//replicate line 0
797  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
798  else
799  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
800  }
801  }
802  //initialize estLineVector
803  if(verbose_opt[0])
804  cout << "initialize estLineVector" << endl;
805  assert(down_opt[0]%2);//window size must be odd
806 
807  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
808  if(iline<0)//replicate line 0
809  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex-1);
810  else
811  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex-1);
812  }
813  statfactory::StatFactory statobs;
814  statobs.setNoDataValues(obsnodata_opt);
815  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
816  //todo: read entire window for uncertReadBuffer...
817  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
818  imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex-1);
819  imgUpdaterUncert.image2geo(0,irow,geox,geoy);
820  if(model_opt.size()==nmodel){
821  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
822  assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
823  imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
824  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
825  }
826  else{
827  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
828  imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
829  }
830  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
831  imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
832 
833  if(modelmask_opt.size()){
834  imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
835  if(modelmask_opt.size()==nmodel)
836  imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
837  else
838  imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
839  }
840 
841  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
842  estLineVector.erase(estLineVector.begin());
843  imgUpdaterEst.readData(estLineBuffer,maxRow,modindex-1);
844  estLineVector.push_back(estLineBuffer);
845  estLineBuffer=estLineVector[down_opt[0]/2];
846 
847  if(update){
848  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
849  obsLineVector.erase(obsLineVector.begin());
850  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
851  obsLineVector.push_back(obsLineBuffer);
852  obsLineBuffer=obsLineVector[down_opt[0]/2];
853 
854  if(observationmask_opt.size())
855  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
856  }
857 
858  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
859  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
860  imgUpdaterEst.image2geo(icol,irow,geox,geoy);
861  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
862  int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
863  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
864  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
865  estWindowBuffer.clear();
866  for(int iline=0;iline<estLineVector.size();++iline){
867  for(int isample=minCol;isample<=maxCol;++isample){
868  assert(isample<estLineVector[iline].size());
869  estWindowBuffer.push_back(estLineVector[iline][isample]);
870  }
871  }
872  if(update){
873  obsWindowBuffer.clear();
874  for(int iline=0;iline<obsLineVector.size();++iline){
875  for(int isample=minCol;isample<=maxCol;++isample){
876  assert(isample<obsLineVector[iline].size());
877  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
878  }
879  }
880  }
881  double estValue=estLineBuffer[icol];
882  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
883  bool model1IsNoData=false;
884  if(modelmask_opt.size())
885  model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
886  lowerCol=modCol-0.5;
887  lowerCol=static_cast<int>(lowerCol);
888  upperCol=modCol+0.5;
889  upperCol=static_cast<int>(upperCol);
890  if(lowerCol<0)
891  lowerCol=0;
892  if(upperCol>=imgReaderModel1.nrOfCol())
893  upperCol=imgReaderModel1.nrOfCol()-1;
894  double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
895  model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
896  if(model_opt.size()==nmodel)
897  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
898  else
899  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
900  bool model2IsNoData=false;
901  if(modelmask_opt.size())
902  model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
903  lowerCol=modCol-0.5;
904  lowerCol=static_cast<int>(lowerCol);
905  upperCol=modCol+0.5;
906  upperCol=static_cast<int>(upperCol);
907  if(lowerCol<0)
908  lowerCol=0;
909  if(upperCol>=imgReaderModel1.nrOfCol())
910  upperCol=imgReaderModel1.nrOfCol()-1;
911  double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
912  model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
913  bool obsIsNoData=false;
914  if(observationmask_opt.size())
915  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
916  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
917 
918  if(imgUpdaterEst.isNoData(estValue)){
919  //we have not found any valid data yet, better here to take the current model value if valid
920  if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
921  estWriteBuffer[icol]=obsnodata_opt[0];
922  uncertWriteBuffer[icol]=uncertNodata_opt[0];
923  //test
924  gainWriteBuffer[icol]=0;
925  }
926  else{
927  estWriteBuffer[icol]=modValue2;
928  uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
929  if(obsmin_opt.size()){
930  if(estWriteBuffer[icol]<obsmin_opt[0])
931  estWriteBuffer[icol]=obsmin_opt[0];
932  }
933  if(obsmax_opt.size()){
934  if(estWriteBuffer[icol]>obsmax_opt[0])
935  estWriteBuffer[icol]=obsmax_opt[0];
936  if(uncertWriteBuffer[icol]>obsmax_opt[0])
937  uncertWriteBuffer[icol]=obsmax_opt[0];
938  }
939  //test
940  gainWriteBuffer[icol]=0;
941  }
942  }
943  else{//previous estimate is valid
944  double estMeanValue=0;
945  double estVarValue=0;
946  statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
947  double nvalid=0;
948  //time update
949  double processNoiseVariance=processNoise_opt[0]*estVarValue;
950  //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
951 
952  if(model1IsNoData||model2IsNoData){
953  estWriteBuffer[icol]=estValue;
954  // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
955  //todo: check following line if makes sense
956  uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
957  }
958  else{//model is good
959  double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
960  estWriteBuffer[icol]=estValue*modRatio;
961  uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
962  }
963  if(obsmin_opt.size()){
964  if(estWriteBuffer[icol]<obsmin_opt[0])
965  estWriteBuffer[icol]=obsmin_opt[0];
966  }
967  if(obsmax_opt.size()){
968  if(estWriteBuffer[icol]>obsmax_opt[0])
969  estWriteBuffer[icol]=obsmax_opt[0];
970  if(uncertWriteBuffer[icol]>obsmax_opt[0])
971  uncertWriteBuffer[icol]=obsmax_opt[0];
972  }
973  }
974  //measurement update
975  if(update&&!obsIsNoData){
976  double kalmanGain=1;
977  if(!model2IsNoData){//model is valid
978  statfactory::StatFactory statobs;
979  statobs.setNoDataValues(obsnodata_opt);
980  double obsMeanValue=0;
981  double obsVarValue=0;
982  double difference=0;
983  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
984  difference=obsMeanValue-modValue2;
985  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
986  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
987 
988  if(errObs<eps_opt[0])
989  errObs=eps_opt[0];
990  double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
991 
992  if(errorCovariance+errObs>eps_opt[0])
993  kalmanGain=errorCovariance/(errorCovariance+errObs);
994  else
995  kalmanGain=1;
996  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
997  uncertWriteBuffer[icol]*=(1-kalmanGain);
998  if(obsmin_opt.size()){
999  if(estWriteBuffer[icol]<obsmin_opt[0])
1000  estWriteBuffer[icol]=obsmin_opt[0];
1001  }
1002  if(obsmax_opt.size()){
1003  if(estWriteBuffer[icol]>obsmax_opt[0])
1004  estWriteBuffer[icol]=obsmax_opt[0];
1005  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1006  uncertWriteBuffer[icol]=obsmax_opt[0];
1007  }
1008  }
1009  if(kalmanGain>=1)
1010  kalmanGain=1;
1011  //test
1012  gainWriteBuffer[icol]=kalmanGain;
1013  }
1014  }
1015  }
1016 
1017  //test
1018  if(gain_opt.size())
1019  imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
1020  imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
1021  imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
1022  progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
1023  pfnProgress(progress,pszMessage,pProgressArg);
1024  }
1025  }
1026 
1027  //must close writers to ensure flush
1028  imgUpdaterEst.close();
1029  imgUpdaterUncert.close();
1030 
1031  if(update){
1032  if(observation_opt.size()==nobs)
1033  imgReaderObs.close();
1034  if(observationmask_opt.size()==nobs)
1035  imgReaderObsMask.close();
1036  ++obsindex;
1037  }
1038  if(model_opt.size()==nmodel){
1039  imgReaderModel1.close();
1040  imgReaderModel2.close();
1041  }
1042  if(modelmask_opt.size()==nmodel){
1043  imgReaderModel1Mask.close();
1044  imgReaderModel2Mask.close();
1045  }
1046  }
1047  //test
1048  if(gain_opt.size())
1049  imgWriterGain.close();
1050  }
1051  }
1052  catch(string errorString){
1053  cerr << errorString << endl;
1054  exit(1);
1055  }
1056  catch(...){
1057  cerr << "Error in forward direction " << endl;
1058  exit(2);
1059  }
1060  try{
1061  if(find(direction_opt.begin(),direction_opt.end(),"backward")!=direction_opt.end()){
1063  cout << "Running backward model" << endl;
1064  obsindex=relobsindex.size()-1;
1065  if(verbose_opt[0])
1066  cout << "Opening image " << outputbw_opt[0] << " for writing " << endl;
1067 
1068  ImgWriterGdal imgWriterEst;
1069  ImgWriterGdal imgWriterUncert;
1070  imgWriterEst.open(outputbw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1071  imgWriterEst.setProjectionProj4(projection_opt[0]);
1072  imgWriterEst.setGeoTransform(geotransform);
1073  imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
1074  imgWriterUncert.open(uncertbw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1075  imgWriterUncert.setProjectionProj4(projection_opt[0]);
1076  imgWriterUncert.setGeoTransform(geotransform);
1077 
1078  try{
1079  // //test
1080  // if(gain_opt.size()){
1081  // imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
1082  // imgWriterGain.setProjectionProj4(projection_opt[0]);
1083  // imgWriterGain.setGeoTransform(geotransform);
1084  // imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
1085  // }
1086 
1087  if(verbose_opt[0]){
1088  cout << "processing time " << tmodel_opt.back() << endl;
1089  if(obsindex<relobsindex.size()){
1090  assert(tmodel_opt.size()>relobsindex[obsindex]);
1091  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
1092  }
1093  else
1094  cout << "There is no next observation" << endl;
1095  }
1096  if(model_opt.size()==nmodel){
1097  imgReaderModel1.open(model_opt.back());
1098  imgReaderModel1.setNoData(modnodata_opt);
1099  }
1100  if(modelmask_opt.size()==nmodel){
1101  imgReaderModel1Mask.open(modelmask_opt[0]);
1102  imgReaderModel1Mask.setNoData(msknodata_opt);
1103  }
1104  }
1105  catch(string errorString){
1106  cerr << errorString << endl;
1107  }
1108  catch(...){
1109  cerr << "Error opening file " << model_opt[0] << endl;
1110  }
1111 
1112  double modRow=0;
1113  double modCol=0;
1114  double lowerCol=0;
1115  double upperCol=0;
1116  RESAMPLE theResample=BILINEAR;
1117 
1118  if(relobsindex.back()<nmodel-1){//initialize output_opt.back() as last model
1119  //write last model as output
1120  if(verbose_opt[0])
1121  cout << "write last model as output" << endl;
1122  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
1123  vector<double> estReadBuffer;
1124  vector<double> lineModelMask;
1125  vector<double> estWriteBuffer(ncol);
1126  vector<double> uncertWriteBuffer(ncol);
1127  // //test
1128  // vector<double> gainWriteBuffer(ncol);
1129  try{
1130  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
1131  imgWriterEst.image2geo(0,irow,geox,geoy);
1132  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1133  if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
1134  cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
1135  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1136  }
1137  // imgReaderModel1.readData(estReadBuffer,modRow,0,theResample);
1138  int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
1139  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
1140  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
1141  if(modelmask_opt.size())
1142  imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
1143  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
1144  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
1145  imgWriterEst.image2geo(icol,irow,geox,geoy);
1146  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1147  if(lineModelMask.size()>modCol){
1148  if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
1149  estWriteBuffer[icol]=obsnodata_opt[0];
1150  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1151  //test
1152  // gainWriteBuffer[icol]=obsnodata_opt[0];
1153  continue;
1154  }
1155  }
1156  lowerCol=modCol-0.5;
1157  lowerCol=static_cast<int>(lowerCol);
1158  upperCol=modCol+0.5;
1159  upperCol=static_cast<int>(upperCol);
1160  if(lowerCol<0)
1161  lowerCol=0;
1162  if(upperCol>=imgReaderModel1.nrOfCol())
1163  upperCol=imgReaderModel1.nrOfCol()-1;
1164  double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
1165  // double modValue=estReadBuffer[modCol];
1166  if(imgReaderModel1.isNoData(modValue)){
1167  estWriteBuffer[icol]=obsnodata_opt[0];
1168  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1169  //test
1170  // gainWriteBuffer[icol]=obsnodata_opt[0];
1171  continue;
1172  }
1173  estWriteBuffer[icol]=modValue;
1174  if(obsmin_opt.size()){
1175  if(estWriteBuffer[icol]<obsmin_opt[0])
1176  estWriteBuffer[icol]=obsmin_opt[0];
1177  }
1178  if(obsmax_opt.size()){
1179  if(estWriteBuffer[icol]>obsmax_opt[0])
1180  estWriteBuffer[icol]=obsmax_opt[0];
1181  }
1182  uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
1183  //test
1184  // gainWriteBuffer[icol]=0;
1185  }
1186  }
1187  imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
1188  imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
1189  // //test
1190  // if(gain_opt.size())
1191  // imgWriterGain.writeData(gainWriteBuffer,irow,0);
1192  }
1193  }
1194  catch(string errorString){
1195  cerr << errorString << endl;
1196  }
1197  catch(...){
1198  cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
1199  }
1200  }
1201  }
1202  else{//we have a measurement at end time
1203  if(verbose_opt[0])
1204  cout << "we have a measurement at end time" << endl;
1205  if(observation_opt.size()==nobs){
1206  imgReaderObs.open(observation_opt.back());
1207  imgReaderObs.setNoData(obsnodata_opt);
1208  }
1209  if(observationmask_opt.size()==nobs){
1210  imgReaderObsMask.open(observationmask_opt.back());
1211  imgReaderObsMask.setNoData(msknodata_opt);
1212  }
1213  imgReaderObs.getGeoTransform(geotransform);
1214 
1215  vector< vector<double> > obsLineVector(down_opt[0]);
1216  vector<double> obsLineBuffer;
1217  vector<double> obsMaskLineBuffer;
1218  vector<double> modelMaskLineBuffer;
1219  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
1220  vector<double> estReadBuffer;
1221  vector<double> estWriteBuffer(ncol);
1222  vector<double> uncertWriteBuffer(ncol);
1223  vector<double> uncertObsLineBuffer;
1224  // //test
1225  // vector<double> gainWriteBuffer(ncol);
1226 
1227  if(verbose_opt[0])
1228  cout << "initialize obsLineVector" << endl;
1229  assert(down_opt[0]%2);//window size must be odd
1230  int readObsBand=(observation_opt.size()==nobs)? 0:nobs-1;
1231  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:nobs-1;
1232  int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
1233  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:nmodel-1;
1234  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
1235  if(iline<0)//replicate line 0
1236  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
1237  else
1238  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
1239  }
1240  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
1241  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
1242  imgWriterEst.image2geo(0,irow,geox,geoy);
1243  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1244  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1245  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
1246  if(modelmask_opt.size())
1247  imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
1248  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
1249  obsLineVector.erase(obsLineVector.begin());
1250  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
1251  obsLineVector.push_back(obsLineBuffer);
1252 
1253  if(observationmask_opt.size())
1254  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
1255 
1256  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
1257  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
1258  imgWriterEst.image2geo(icol,irow,geox,geoy);
1259  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1260  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1261  bool modelIsNoData=false;
1262  if(modelmask_opt.size())
1263  modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
1264  lowerCol=modCol-0.5;
1265  lowerCol=static_cast<int>(lowerCol);
1266  upperCol=modCol+0.5;
1267  upperCol=static_cast<int>(upperCol);
1268  if(lowerCol<0)
1269  lowerCol=0;
1270  if(upperCol>=imgReaderModel1.nrOfCol())
1271  upperCol=imgReaderModel1.nrOfCol()-1;
1272  double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
1273  // double modValue=estReadBuffer[modCol];
1274  double errMod=uncertModel_opt[0];//*stdDev*stdDev;
1275  modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
1276  bool obsIsNoData=false;
1277  if(observationmask_opt.size())
1278  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
1279  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
1280  if(modelIsNoData){//model is nodata: retain observation
1281  if(obsIsNoData){//both model and observation nodata
1282  estWriteBuffer[icol]=obsnodata_opt[0];
1283  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1284  //test
1285  // gainWriteBuffer[icol]=obsnodata_opt[0];
1286  }
1287  else{
1288  estWriteBuffer[icol]=obsLineBuffer[icol];
1289  if(obsmin_opt.size()){
1290  if(estWriteBuffer[icol]<obsmin_opt[0])
1291  estWriteBuffer[icol]=obsmin_opt[0];
1292  }
1293  if(obsmax_opt.size()){
1294  if(estWriteBuffer[icol]>obsmax_opt[0])
1295  estWriteBuffer[icol]=obsmax_opt[0];
1296  }
1297  uncertWriteBuffer[icol]=uncertObs_opt[0];
1298  }
1299  }
1300  else{//model is valid: calculate estimate from model
1301  estWriteBuffer[icol]=modValue;
1302  uncertWriteBuffer[icol]=errMod;//in case observation is not valid
1303  //test
1304  // gainWriteBuffer[icol]=0;
1305  }
1306  //measurement update
1307  if(!obsIsNoData){
1308  // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
1309  double kalmanGain=1;
1310  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
1311  int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
1312  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
1313  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
1314  obsWindowBuffer.clear();
1315  for(int iline=0;iline<obsLineVector.size();++iline){
1316  for(int isample=minCol;isample<=maxCol;++isample){
1317  assert(isample<obsLineVector[iline].size());
1318  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
1319  }
1320  }
1321  if(!modelIsNoData){//model is valid
1322  statfactory::StatFactory statobs;
1323  statobs.setNoDataValues(obsnodata_opt);
1324  double obsMeanValue=0;
1325  double obsVarValue=0;
1326  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
1327  double difference=0;
1328  difference=obsMeanValue-modValue;
1329  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
1330  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
1331  // double errorCovariance=errMod;
1332  double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
1333  if(errorCovariance+errObs>eps_opt[0])
1334  kalmanGain=errorCovariance/(errorCovariance+errObs);
1335  else
1336  kalmanGain=1;
1337  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
1338  if(obsmin_opt.size()){
1339  if(estWriteBuffer[icol]<obsmin_opt[0])
1340  estWriteBuffer[icol]=obsmin_opt[0];
1341  }
1342  if(obsmax_opt.size()){
1343  if(estWriteBuffer[icol]>obsmax_opt[0])
1344  estWriteBuffer[icol]=obsmax_opt[0];
1345  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1346  uncertWriteBuffer[icol]=obsmax_opt[0];
1347  }
1348  }
1349  if(kalmanGain>=1)
1350  kalmanGain=1;
1351  //test
1352  // gainWriteBuffer[icol]=kalmanGain;
1353  }
1354  }
1355  }
1356  imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
1357  imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
1358  // //test
1359  // if(gain_opt.size())
1360  // imgWriterGain.writeData(gainWriteBuffer,irow,0);
1361  }
1362  }
1363  if(observation_opt.size()==nobs)
1364  imgReaderObs.close();
1365  if(observationmask_opt.size()==nobs)
1366  imgReaderObsMask.close();
1367  --obsindex;
1368  }
1369 
1370  if(model_opt.size()==nmodel)
1371  imgReaderModel1.close();
1372  if(modelmask_opt.size()==nmodel)
1373  imgReaderModel1Mask.close();
1374  imgWriterEst.close();
1375  imgWriterUncert.close();
1376 
1377  ImgUpdaterGdal imgUpdaterEst;
1378  ImgUpdaterGdal imgUpdaterUncert;
1379  for(int modindex=nmodel-2;modindex>=0;--modindex){
1380  imgUpdaterEst.open(outputbw_opt[0]);
1381  imgUpdaterEst.setNoData(obsnodata_opt);
1382  imgUpdaterUncert.open(uncertbw_opt[0]);
1383  if(verbose_opt[0]){
1384  cout << "processing time " << tmodel_opt[modindex] << endl;
1385  if(obsindex<relobsindex.size())
1386  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
1387  else
1388  cout << "There is no next observation" << endl;
1389  }
1390 
1391  //calculate regression between two subsequence model inputs
1392  if(model_opt.size()==nmodel){
1393  imgReaderModel1.open(model_opt[modindex+1]);
1394  imgReaderModel1.setNoData(modnodata_opt);
1395  imgReaderModel2.open(model_opt[modindex]);
1396  imgReaderModel2.setNoData(modnodata_opt);
1397  }
1398  if(modelmask_opt.size()==nmodel){
1399  imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
1400  imgReaderModel1Mask.setNoData(msknodata_opt);
1401  imgReaderModel2Mask.open(modelmask_opt[modindex]);
1402  imgReaderModel2Mask.setNoData(msknodata_opt);
1403  }
1404 
1405  pfnProgress(progress,pszMessage,pProgressArg);
1406 
1407  bool update=false;
1408  if(obsindex<relobsindex.size()){
1409  update=(relobsindex[obsindex]==modindex);
1410  }
1411  if(update){
1412  if(observation_opt.size()==nobs){
1413  if(verbose_opt[0])
1414  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
1415  imgReaderObs.open(observation_opt[obsindex]);
1416  imgReaderObs.getGeoTransform(geotransform);
1417  imgReaderObs.setNoData(obsnodata_opt);
1418  }
1419  if(observationmask_opt.size()==nobs){
1420  imgReaderObsMask.open(observationmask_opt[obsindex]);
1421  imgReaderObsMask.setNoData(msknodata_opt);
1422  }
1423  }
1424  //prediction (also to fill cloudy pixels in update mode)
1425  string input;
1426  input=outputbw_opt[0];
1427 
1428  vector< vector<double> > obsLineVector(down_opt[0]);
1429  vector<double> obsLineBuffer;
1430  vector<double> obsMaskLineBuffer;
1431  vector<double> model1MaskLineBuffer;
1432  vector<double> model2MaskLineBuffer;
1433  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
1434  vector<double> model1LineBuffer;
1435  vector<double> model2LineBuffer;
1436  vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
1437  vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
1438  vector<double> uncertObsLineBuffer;
1439  vector< vector<double> > estLineVector(down_opt[0]);
1440  vector<double> estLineBuffer;
1441  vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
1442  vector<double> uncertReadBuffer;
1443  vector<double> estWriteBuffer(ncol);
1444  vector<double> uncertWriteBuffer(ncol);
1445  //test
1446  // vector<double> gainWriteBuffer(ncol);
1447 
1448  int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
1449  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
1450  int readModel1Band=(model_opt.size()==nmodel)? 0:modindex+1;
1451  int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
1452  int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex+1;
1453  int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
1454 
1455  //initialize obsLineVector
1456  if(update){
1457  if(verbose_opt[0])
1458  cout << "initialize obsLineVector" << endl;
1459  assert(down_opt[0]%2);//window size must be odd
1460  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
1461  if(iline<0)//replicate line 0
1462  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
1463  else
1464  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
1465  }
1466  }
1467  //initialize estLineVector
1468  if(verbose_opt[0])
1469  cout << "initialize estLineVector" << endl;
1470  assert(down_opt[0]%2);//window size must be odd
1471 
1472  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
1473  if(iline<0)//replicate line 0
1474  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex+1);
1475  else
1476  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex+1);
1477  }
1478  statfactory::StatFactory statobs;
1479  statobs.setNoDataValues(obsnodata_opt);
1480 
1481  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
1482  //todo: read entire window for uncertReadBuffer...
1483  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
1484  imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex+1);
1485  imgUpdaterUncert.image2geo(0,irow,geox,geoy);
1486  if(model_opt.size()==nmodel){
1487  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
1488  assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
1489  imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
1490  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1491  }
1492  else{
1493  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1494  imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
1495  }
1496 
1497  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1498  imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
1499  if(modelmask_opt.size()){
1500  imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
1501  if(modelmask_opt.size()==nmodel)
1502  imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
1503  else
1504  imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
1505  }
1506  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
1507  estLineVector.erase(estLineVector.begin());
1508  imgUpdaterEst.readData(estLineBuffer,maxRow,modindex+1);
1509  estLineVector.push_back(estLineBuffer);
1510  estLineBuffer=estLineVector[down_opt[0]/2];
1511 
1512  if(update){
1513  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
1514  obsLineVector.erase(obsLineVector.begin());
1515  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
1516  obsLineVector.push_back(obsLineBuffer);
1517  obsLineBuffer=obsLineVector[down_opt[0]/2];
1518 
1519  if(observationmask_opt.size())
1520  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsBand);
1521  }
1522  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
1523  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
1524  imgUpdaterEst.image2geo(icol,irow,geox,geoy);
1525  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
1526  int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
1527  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
1528  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
1529  estWindowBuffer.clear();
1530  for(int iline=0;iline<estLineVector.size();++iline){
1531  for(int isample=minCol;isample<=maxCol;++isample){
1532  assert(isample<estLineVector[iline].size());
1533  estWindowBuffer.push_back(estLineVector[iline][isample]);
1534  }
1535  }
1536  if(update){
1537  obsWindowBuffer.clear();
1538  for(int iline=0;iline<obsLineVector.size();++iline){
1539  for(int isample=minCol;isample<=maxCol;++isample){
1540  assert(isample<obsLineVector[iline].size());
1541  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
1542  }
1543  }
1544  }
1545 
1546  double estValue=estLineBuffer[icol];
1547  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1548  bool model1IsNoData=false;
1549 
1550  if(modelmask_opt.size())
1551  model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
1552 
1553  lowerCol=modCol-0.5;
1554  lowerCol=static_cast<int>(lowerCol);
1555  upperCol=modCol+0.5;
1556  upperCol=static_cast<int>(upperCol);
1557  if(lowerCol<0)
1558  lowerCol=0;
1559  if(upperCol>=imgReaderModel1.nrOfCol())
1560  upperCol=imgReaderModel1.nrOfCol()-1;
1561  double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
1562  model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
1563  if(model_opt.size()==nmodel)
1564  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
1565  else
1566  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1567  bool model2IsNoData=false;
1568 
1569  if(modelmask_opt.size())
1570  model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
1571  lowerCol=modCol-0.5;
1572  lowerCol=static_cast<int>(lowerCol);
1573  upperCol=modCol+0.5;
1574  upperCol=static_cast<int>(upperCol);
1575  if(lowerCol<0)
1576  lowerCol=0;
1577  if(upperCol>=imgReaderModel1.nrOfCol())
1578  upperCol=imgReaderModel1.nrOfCol()-1;
1579  double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
1580  model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
1581  bool obsIsNoData=false;
1582  if(observationmask_opt.size())
1583  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
1584  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
1585 
1586  if(imgUpdaterEst.isNoData(estValue)){
1587  //we have not found any valid data yet, better here to take the current model value if valid
1588  if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
1589  estWriteBuffer[icol]=obsnodata_opt[0];
1590  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1591  //test
1592  // gainWriteBuffer[icol]=0;
1593  }
1594  else{
1595  estWriteBuffer[icol]=modValue2;
1596  uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
1597  if(obsmin_opt.size()){
1598  if(estWriteBuffer[icol]<obsmin_opt[0])
1599  estWriteBuffer[icol]=obsmin_opt[0];
1600  }
1601  if(obsmax_opt.size()){
1602  if(estWriteBuffer[icol]>obsmax_opt[0])
1603  estWriteBuffer[icol]=obsmax_opt[0];
1604  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1605  uncertWriteBuffer[icol]=obsmax_opt[0];
1606  }
1607  //test
1608  // gainWriteBuffer[icol]=0;
1609  }
1610  }
1611  else{//previous estimate is valid
1612  double estMeanValue=0;
1613  double estVarValue=0;
1614  statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
1615  double nvalid=0;
1616  //time update
1617  double processNoiseVariance=processNoise_opt[0]*estVarValue;
1618  //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
1619 
1620  if(model1IsNoData||model2IsNoData){
1621  estWriteBuffer[icol]=estValue;
1622  // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
1623  //todo: check following line if makes sense
1624  uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
1625  }
1626  else{//model is good
1627  double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
1628  estWriteBuffer[icol]=estValue*modRatio;
1629  uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
1630  }
1631  if(obsmin_opt.size()){
1632  if(estWriteBuffer[icol]<obsmin_opt[0])
1633  estWriteBuffer[icol]=obsmin_opt[0];
1634  }
1635  if(obsmax_opt.size()){
1636  if(estWriteBuffer[icol]>obsmax_opt[0])
1637  estWriteBuffer[icol]=obsmax_opt[0];
1638  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1639  uncertWriteBuffer[icol]=obsmax_opt[0];
1640  }
1641  }
1642  //measurement update
1643  if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
1644  double kalmanGain=1;
1645  if(!imgReaderModel1.isNoData(modValue2)){//model is valid
1646  statfactory::StatFactory statobs;
1647  statobs.setNoDataValues(obsnodata_opt);
1648  double obsMeanValue=0;
1649  double obsVarValue=0;
1650  double difference=0;
1651  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
1652  difference=obsMeanValue-modValue2;
1653  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
1654  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
1655 
1656  if(errObs<eps_opt[0])
1657  errObs=eps_opt[0];
1658  double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
1659 
1660  if(errorCovariance+errObs>eps_opt[0])
1661  kalmanGain=errorCovariance/(errorCovariance+errObs);
1662  else
1663  kalmanGain=1;
1664  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
1665  uncertWriteBuffer[icol]*=(1-kalmanGain);
1666  if(obsmin_opt.size()){
1667  if(estWriteBuffer[icol]<obsmin_opt[0])
1668  estWriteBuffer[icol]=obsmin_opt[0];
1669  }
1670  if(obsmax_opt.size()){
1671  if(estWriteBuffer[icol]>obsmax_opt[0])
1672  estWriteBuffer[icol]=obsmax_opt[0];
1673  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1674  uncertWriteBuffer[icol]=obsmax_opt[0];
1675  }
1676  }
1677  if(kalmanGain>=1)
1678  kalmanGain=1;
1679  //test
1680  // gainWriteBuffer[icol]=kalmanGain;
1681  }
1682  }
1683  }
1684  // //test
1685  // if(gain_opt.size())
1686  // imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
1687  imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
1688  imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
1689  progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
1690  pfnProgress(progress,pszMessage,pProgressArg);
1691  }
1692  }
1693  //must close writers to ensure flush
1694  imgUpdaterEst.close();
1695  imgUpdaterUncert.close();
1696 
1697  if(update){
1698  if(observation_opt.size()==nobs)
1699  imgReaderObs.close();
1700  if(observationmask_opt.size()==nobs)
1701  imgReaderObsMask.close();
1702  --obsindex;
1703  }
1704  if(model_opt.size()==nmodel){
1705  imgReaderModel1.close();
1706  imgReaderModel2.close();
1707  }
1708  if(modelmask_opt.size()==nmodel){
1709  imgReaderModel1Mask.close();
1710  imgReaderModel2Mask.close();
1711  }
1712  }
1713  // //test
1714  // if(gain_opt.size())
1715  // imgWriterGain.close();
1716  }
1717  }
1718  catch(string errorString){
1719  cerr << errorString << endl;
1720  exit(1);
1721  }
1722  catch(...){
1723  cerr << "Error in backward direction " << endl;
1724  exit(2);
1725  }
1726  if(find(direction_opt.begin(),direction_opt.end(),"smooth")!=direction_opt.end()){
1728  cout << "Running smooth model" << endl;
1729  obsindex=0;
1730 
1731  ImgReaderGdal imgReaderForward(outputfw_opt[0]);
1732  ImgReaderGdal imgReaderBackward(outputbw_opt[0]);
1733  ImgReaderGdal imgReaderForwardUncert(uncertfw_opt[0]);
1734  ImgReaderGdal imgReaderBackwardUncert(uncertbw_opt[0]);
1735  imgReaderForward.setNoData(obsnodata_opt);
1736  imgReaderBackward.setNoData(obsnodata_opt);
1737 
1738  assert(imgReaderForward.nrOfBand()==nmodel);
1739  assert(imgReaderForwardUncert.nrOfBand()==nmodel);
1740  assert(imgReaderBackward.nrOfBand()==nmodel);
1741  assert(imgReaderBackwardUncert.nrOfBand()==nmodel);
1742  ImgWriterGdal imgWriterEst;
1743  imgWriterEst.setNoData(obsnodata_opt);
1744  ImgWriterGdal imgWriterUncert;
1745  imgWriterEst.open(outputfb_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1746  imgWriterEst.setProjectionProj4(projection_opt[0]);
1747  imgWriterEst.setGeoTransform(geotransform);
1748  imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
1749 
1750  imgWriterUncert.open(uncertfb_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1751  imgWriterUncert.setProjectionProj4(projection_opt[0]);
1752  imgWriterUncert.setGeoTransform(geotransform);
1753  for(int modindex=0;modindex<nmodel;++modindex){
1754  if(verbose_opt[0]){
1755  cout << "processing time " << tmodel_opt[modindex] << endl;
1756  if(obsindex<relobsindex.size())
1757  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
1758  else
1759  cout << "There is no next observation" << endl;
1760  }
1761  vector<double> estForwardBuffer;
1762  vector<double> estBackwardBuffer;
1763  vector<double> uncertObsLineBuffer;
1764  vector<double> uncertForwardBuffer;
1765  vector<double> uncertBackwardBuffer;
1766  vector<double> uncertReadBuffer;
1767  vector<double> estWriteBuffer(ncol);
1768  vector<double> uncertWriteBuffer(ncol);
1769 
1770  bool update=false;
1771  if(obsindex<relobsindex.size()){
1772  update=(relobsindex[obsindex]==modindex);
1773  }
1774 
1775  int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
1776 
1777  if(update){
1778  if(observation_opt.size()==nobs){
1779  if(verbose_opt[0])
1780  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
1781  imgReaderObs.open(observation_opt[obsindex]);
1782  imgReaderObs.setNoData(obsnodata_opt);
1783  imgReaderObs.getGeoTransform(geotransform);
1784  }
1785  if(observationmask_opt.size()==nobs){
1786  imgReaderObsMask.open(observationmask_opt[obsindex]);
1787  imgReaderObsMask.setNoData(msknodata_opt);
1788  }
1789  }
1790 
1791  pfnProgress(progress,pszMessage,pProgressArg);
1792 
1793  for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
1794  assert(irow<imgReaderForward.nrOfRow());
1795  assert(irow<imgReaderBackward.nrOfRow());
1796  imgReaderForward.readData(estForwardBuffer,irow,modindex);
1797  imgReaderBackward.readData(estBackwardBuffer,irow,modindex);
1798  imgReaderForwardUncert.readData(uncertForwardBuffer,irow,modindex);
1799  imgReaderBackwardUncert.readData(uncertBackwardBuffer,irow,modindex);
1800 
1801  if(update){
1802  if(observation_opt.size()==nobs)
1803  imgReaderObs.readData(estWriteBuffer,irow,readObsBand);
1804  if(observationmask_opt.size())
1805  imgReaderObsMask.readData(uncertObsLineBuffer,irow,readObsBand);
1806  }
1807 
1808  for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
1809  imgWriterEst.image2geo(icol,irow,geox,geoy);
1810  double A=estForwardBuffer[icol];
1811  double B=estBackwardBuffer[icol];
1812  double C=uncertForwardBuffer[icol];
1813  double D=uncertBackwardBuffer[icol];
1814  double uncertObs=uncertObs_opt[0];
1815 
1816  double noemer=(C+D);
1817  //todo: consistently check for division by zero...
1818  if(imgReaderForward.isNoData(A)&&imgReaderBackward.isNoData(B)){
1819  estWriteBuffer[icol]=obsnodata_opt[0];
1820  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1821  }
1822  else if(imgReaderForward.isNoData(A)){
1823  estWriteBuffer[icol]=B;
1824  uncertWriteBuffer[icol]=uncertBackwardBuffer[icol];
1825  }
1826  else if(imgReaderForward.isNoData(B)){
1827  estWriteBuffer[icol]=A;
1828  uncertWriteBuffer[icol]=uncertForwardBuffer[icol];
1829  }
1830  else{
1831  if(noemer<eps_opt[0]){//simple average if both uncertainties are ~>0
1832  estWriteBuffer[icol]=0.5*(A+B);
1833  uncertWriteBuffer[icol]=eps_opt[0];
1834  }
1835  else{
1836  estWriteBuffer[icol]=(A*D+B*C)/noemer;
1837  uncertWriteBuffer[icol]=C*D/noemer;
1838  if(obsmin_opt.size()){
1839  if(estWriteBuffer[icol]<obsmin_opt[0])
1840  estWriteBuffer[icol]=obsmin_opt[0];
1841  }
1842  if(obsmax_opt.size()){
1843  if(estWriteBuffer[icol]>obsmax_opt[0])
1844  estWriteBuffer[icol]=obsmax_opt[0];
1845  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1846  uncertWriteBuffer[icol]=obsmax_opt[0];
1847  }
1848  }
1849  }
1850  }
1851  imgWriterEst.writeData(estWriteBuffer,irow,modindex);
1852  imgWriterUncert.writeData(uncertWriteBuffer,irow,modindex);
1853  progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
1854  pfnProgress(progress,pszMessage,pProgressArg);
1855  }
1856  if(update){
1857  if(observation_opt.size()==nobs)
1858  imgReaderObs.close();
1859  ++obsindex;
1860  }
1861  }
1862  imgReaderForward.close();
1863  imgReaderBackward.close();
1864  imgWriterEst.close();
1865  imgWriterUncert.close();
1866  }
1867  if(observation_opt.size()<nobs)
1868  imgReaderObs.close();
1869  if(model_opt.size()<nmodel)
1870  imgReaderModel1.close();
1871 }
1872 
void close(void)
Close the image.
std::string getFileName() const
Get the filename of this dataset.
Definition: ImgRasterGdal.h:96
bool isNoData(double value) const
Check if value is nodata in this dataset.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
CPLErr setProjectionProj4(const std::string &projection)
Set the projection for this dataset from user input (supports epsg: format) ...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y) ...
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
int setNoData(const std::vector< double > nodata)
Set the no data values of this dataset using a standard template library (stl) vector as input...
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
double getDeltaX(void) const
Get the pixel cell spacing in x.
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98