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  double down_f=resModel/resObs;
328  // round up to the next integer value, unless we already have
329  // an integer (or a value close enough)
330  int down=static_cast<int>(ceil(down_f - eps_opt[0]));
331  // the downsampling factor should be odd
332 
333  if(!(down%2))
334  down+=1;
335  down_opt.push_back(down);
336  }
337 
338  int obsindex=0;
339 
340  const char* pszMessage;
341  void* pProgressArg=NULL;
342  GDALProgressFunc pfnProgress=GDALTermProgress;
343  double progress=0;
344 
345  double errObs=uncertNodata_opt[0];//start with high initial value in case we do not have first observation at time 0
346 
347  while(tmodel_opt.size()<nmodel)
348  tmodel_opt.push_back(tmodel_opt.size());
349  try{
350  if(tobservation_opt.size()<nobs){
351  if(nobs==nmodel){
352  while(tobservation_opt.size()<nobs)
353  tobservation_opt.push_back(tobservation_opt.size());
354  }
355  else{
356  ostringstream errorStream;
357  errorStream << "Error: please provide time sequence for observation using option tobs" << endl;
358  throw(errorStream.str());
359  }
360  }
361  }
362  catch(string errorString){
363  std::cout << errorString << std::endl;
364  exit(1);
365  }
366 
367  vector<int> relobsindex;
368 
369  for(int tindex=0;tindex<tobservation_opt.size();++tindex){
370  vector<int>::iterator modit;
371  modit=upper_bound(tmodel_opt.begin(),tmodel_opt.end(),tobservation_opt[tindex]);
372  int relpos=modit-tmodel_opt.begin()-1;
373  assert(relpos>=0);//todo: for now, we assume model is available at time before first measurement
374  relobsindex.push_back(relpos);
375  if(verbose_opt[0]){
376  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();
377  if(observation_opt.size()>tindex)
378  cout << ", filename observation: " << observation_opt[tindex];
379  else
380  cout << ", observation band index: " << tindex;
381  if(model_opt.size()>relpos)
382  cout << ", filename of corresponding model: " << model_opt[relpos] << endl;
383  else
384  cout << ", band index of corresponding model: " << relpos;
385  }
386  }
387 
388  int ndigit=log(1.0*tmodel_opt.back())/log(10.0)+1;
389 
390  double geox=0;
391  double geoy=0;
392 
393  if(model_opt.size()==nmodel)
394  imgReaderModel1.close();
395  if(modelmask_opt.size()==nmodel)
396  imgReaderModel1Mask.close();
397  if(observation_opt.size()==nobs)
398  imgReaderObs.close();
399  if(observationmask_opt.size()==nobs)
400  imgReaderObsMask.close();
401 
402  try{
403  if(find(direction_opt.begin(),direction_opt.end(),"forward")!=direction_opt.end()){
405  cout << "Running forward model" << endl;
406  obsindex=0;
407  if(verbose_opt[0])
408  cout << "Opening image " << outputfw_opt[0] << " for writing " << endl << flush;
409 
410  ImgWriterGdal imgWriterEst;
411  ImgWriterGdal imgWriterUncert;
412  imgWriterEst.open(outputfw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
413  imgWriterEst.setProjectionProj4(projection_opt[0]);
414  imgWriterEst.setGeoTransform(geotransform);
415  imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
416  imgWriterUncert.open(uncertfw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
417  imgWriterUncert.setProjectionProj4(projection_opt[0]);
418  imgWriterUncert.setGeoTransform(geotransform);
419 
420  try{
421  //test
422  if(gain_opt.size()){
423  imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
424  imgWriterGain.setProjectionProj4(projection_opt[0]);
425  imgWriterGain.setGeoTransform(geotransform);
426  imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
427  }
428 
429  if(verbose_opt[0]){
430  cout << "processing time " << tmodel_opt[0] << endl;
431  if(obsindex<relobsindex.size()){
432  assert(tmodel_opt.size()>relobsindex[obsindex]);
433  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
434  }
435  else
436  cout << "There is no next observation" << endl;
437  }
438  if(model_opt.size()==nmodel){
439  imgReaderModel1.open(model_opt[0]);
440  imgReaderModel1.setNoData(modnodata_opt);
441  }
442  if(modelmask_opt.size()==nmodel){
443  imgReaderModel1Mask.open(modelmask_opt[0]);
444  imgReaderModel1Mask.setNoData(msknodata_opt);
445  }
446  }
447  catch(string errorString){
448  cerr << errorString << endl;
449  }
450  catch(...){
451  cerr << "Error opening file " << model_opt[0] << endl;
452  }
453 
454  double modEps=0.00001;
455  double modRow=0;
456  double modCol=0;
457  double lowerCol=0;
458  double upperCol=0;
459  RESAMPLE theResample=BILINEAR;
460 
461  if(relobsindex[0]>0){//initialize output_opt[0] as model[0]
462  //write first model as output
463  if(verbose_opt[0])
464  cout << "write first model as output" << endl;
465  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
466  vector<double> estReadBuffer;
467  vector<double> lineModelMask;
468  vector<double> estWriteBuffer(ncol);
469  vector<double> uncertWriteBuffer(ncol);
470  //test
471  vector<double> gainWriteBuffer(ncol);
472  try{
473  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
474  imgWriterEst.image2geo(0,irow,geox,geoy);
475  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
476  if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
477  cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
478  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
479  }
480 
481  int readModelBand=(model_opt.size()==nmodel)? 0:0;
482  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
483  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
484  if(modelmask_opt.size())
485  imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
486  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
487  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
488  imgWriterEst.image2geo(icol,irow,geox,geoy);
489  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
490  if(modelmask_opt.size()){
491  if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
492  estWriteBuffer[icol]=obsnodata_opt[0];
493  uncertWriteBuffer[icol]=uncertNodata_opt[0];
494  //test
495  gainWriteBuffer[icol]=obsnodata_opt[0];
496  continue;
497  }
498  }
499  // lowerCol=modCol-0.5;
500  // lowerCol=static_cast<int>(lowerCol);
501  // upperCol=modCol+0.5;
502  // upperCol=static_cast<int>(upperCol);
503  lowerCol=static_cast<int>(modCol-0.5+modEps);
504  if(lowerCol<0)
505  lowerCol=0;
506  if(lowerCol>=imgReaderModel1.nrOfCol())
507  lowerCol=imgReaderModel1.nrOfCol()-1;
508  upperCol=lowerCol+1.0;
509  if(upperCol>=imgReaderModel1.nrOfCol())
510  upperCol=imgReaderModel1.nrOfCol()-1;
511  // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
512  double modValue;
513  if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
514  if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
515  modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
516  }
517  else{
518  modValue=estReadBuffer[lowerCol];
519  }
520  }
521  else{
522  modValue=estReadBuffer[upperCol];
523  }
524  // double modValue=estReadBuffer[modCol];
525 
526  if(imgReaderModel1.isNoData(modValue)){
527  estWriteBuffer[icol]=obsnodata_opt[0];
528  uncertWriteBuffer[icol]=uncertNodata_opt[0];
529  //test
530  gainWriteBuffer[icol]=obsnodata_opt[0];
531  continue;
532  }
533  estWriteBuffer[icol]=modValue;
534  if(obsmin_opt.size()){
535  if(estWriteBuffer[icol]<obsmin_opt[0])
536  estWriteBuffer[icol]=obsmin_opt[0];
537  }
538  if(obsmax_opt.size()){
539  if(estWriteBuffer[icol]>obsmax_opt[0])
540  estWriteBuffer[icol]=obsmax_opt[0];
541  }
542  uncertWriteBuffer[icol]=uncertModel_opt[0];
543  //test
544  gainWriteBuffer[icol]=0;
545  }
546  }
547  imgWriterEst.writeData(estWriteBuffer,irow,0);
548  imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
549  //test
550  if(gain_opt.size())
551  imgWriterGain.writeData(gainWriteBuffer,irow,0);
552  }
553  }
554  catch(string errorString){
555  cerr << errorString << endl;
556  }
557  catch(...){
558  cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
559  }
560  }
561  }
562  else{//we have a measurement
563  if(verbose_opt[0])
564  cout << "we have a measurement at initial time" << endl;
565  if(observation_opt.size()==nobs){
566  imgReaderObs.open(observation_opt[0]);
567  imgReaderObs.setNoData(obsnodata_opt);
568  }
569  if(observationmask_opt.size()==nobs){
570  imgReaderObsMask.open(observationmask_opt[0]);
571  imgReaderObsMask.setNoData(msknodata_opt);
572  }
573  imgReaderObs.getGeoTransform(geotransform);
574 
575  vector< vector<double> > obsLineVector(down_opt[0]);
576  vector<double> obsLineBuffer;
577  vector<double> obsMaskLineBuffer;
578  vector<double> modelMaskLineBuffer;
579  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
580  vector<double> estReadBuffer;
581  vector<double> estWriteBuffer(ncol);
582  vector<double> uncertWriteBuffer(ncol);
583  vector<double> uncertObsLineBuffer;
584  //test
585  vector<double> gainWriteBuffer(ncol);
586 
587  if(verbose_opt[0])
588  cout << "initialize obsLineVector" << endl;
589  assert(down_opt[0]%2);//window size must be odd
590  int readObsBand=(observation_opt.size()==nobs)? 0:0;
591  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:0;
592  int readModelBand=(model_opt.size()==nmodel)? 0:0;
593  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
594  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
595  if(iline<0)//replicate line 0
596  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
597  else
598  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
599  }
600  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
601  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
602  imgWriterEst.image2geo(0,irow,geox,geoy);
603  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
604  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
605  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
606  if(modelmask_opt.size())
607  imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
608  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
609  obsLineVector.erase(obsLineVector.begin());
610  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
611  obsLineVector.push_back(obsLineBuffer);
612 
613  if(observationmask_opt.size())
614  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
615 
616  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
617  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
618  imgWriterEst.image2geo(icol,irow,geox,geoy);
619  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
620  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
621  bool modelIsNoData=false;
622  if(modelmask_opt.size())
623  modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
624  // lowerCol=modCol-0.5;
625  // lowerCol=static_cast<int>(lowerCol);
626  // upperCol=modCol+0.5;
627  // upperCol=static_cast<int>(upperCol);
628  lowerCol=static_cast<int>(modCol-0.5+modEps);
629  if(lowerCol<0)
630  lowerCol=0;
631  if(lowerCol>=imgReaderModel1.nrOfCol())
632  lowerCol=imgReaderModel1.nrOfCol()-1;
633  upperCol=lowerCol+1.0;
634  if(upperCol>=imgReaderModel1.nrOfCol())
635  upperCol=imgReaderModel1.nrOfCol()-1;
636  // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
637  double modValue;
638  if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
639  if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
640  modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
641  }
642  else{
643  modValue=estReadBuffer[lowerCol];
644  }
645  }
646  else{
647  modValue=estReadBuffer[upperCol];
648  }
649  // double modValue=estReadBuffer[modCol];
650 
651  double errMod=uncertModel_opt[0];//*stdDev*stdDev;
652  modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
653  bool obsIsNoData=false;
654  if(observationmask_opt.size())
655  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
656  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
657  if(modelIsNoData){//model is nodata: retain observation
658  if(obsIsNoData){//both model and observation nodata
659  estWriteBuffer[icol]=obsnodata_opt[0];
660  uncertWriteBuffer[icol]=uncertNodata_opt[0];
661  //test
662  gainWriteBuffer[icol]=obsnodata_opt[0];
663  }
664  else{
665  estWriteBuffer[icol]=obsLineBuffer[icol];
666  if(obsmin_opt.size()){
667  if(estWriteBuffer[icol]<obsmin_opt[0])
668  estWriteBuffer[icol]=obsmin_opt[0];
669  }
670  if(obsmax_opt.size()){
671  if(estWriteBuffer[icol]>obsmax_opt[0])
672  estWriteBuffer[icol]=obsmax_opt[0];
673  }
674  uncertWriteBuffer[icol]=uncertObs_opt[0];
675  }
676  }
677  else{//model is valid: calculate estimate from model
678  estWriteBuffer[icol]=modValue;
679  uncertWriteBuffer[icol]=errMod;//in case observation is not valid
680  //test
681  gainWriteBuffer[icol]=0;
682  }
683  //measurement update
684  if(!obsIsNoData){
685  // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
686  double kalmanGain=1;
687  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
688  int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
689  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
690  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
691  obsWindowBuffer.clear();
692  for(int iline=0;iline<obsLineVector.size();++iline){
693  for(int isample=minCol;isample<=maxCol;++isample){
694  assert(isample<obsLineVector[iline].size());
695  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
696  }
697  }
698  if(!modelIsNoData){//model is valid
699  statfactory::StatFactory statobs;
700  statobs.setNoDataValues(obsnodata_opt);
701  double obsMeanValue=0;
702  double obsVarValue=0;
703  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
704  double difference=0;
705  difference=obsMeanValue-modValue;
706  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
707  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
708  // double errorCovariance=errMod;
709  double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
710  if(errorCovariance+errObs>eps_opt[0])
711  kalmanGain=errorCovariance/(errorCovariance+errObs);
712  else
713  kalmanGain=1;
714  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
715  if(obsmin_opt.size()){
716  if(estWriteBuffer[icol]<obsmin_opt[0])
717  estWriteBuffer[icol]=obsmin_opt[0];
718  }
719  if(obsmax_opt.size()){
720  if(estWriteBuffer[icol]>obsmax_opt[0])
721  estWriteBuffer[icol]=obsmax_opt[0];
722  if(uncertWriteBuffer[icol]>obsmax_opt[0])
723  uncertWriteBuffer[icol]=obsmax_opt[0];
724  }
725  }
726  if(kalmanGain>=1)
727  kalmanGain=1;
728  //test
729  gainWriteBuffer[icol]=kalmanGain;
730  }
731  }
732  }
733  imgWriterEst.writeData(estWriteBuffer,irow,0);
734  imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
735  //test
736  if(gain_opt.size())
737  imgWriterGain.writeData(gainWriteBuffer,irow,0);
738  }
739  }
740  if(observation_opt.size()==nobs)
741  imgReaderObs.close();
742  if(observationmask_opt.size()==nobs)
743  imgReaderObsMask.close();
744  ++obsindex;
745  }
746  if(model_opt.size()==nmodel)
747  imgReaderModel1.close();
748  if(modelmask_opt.size()==nmodel)
749  imgReaderModel1Mask.close();
750  imgWriterEst.close();
751  imgWriterUncert.close();
752 
753  ImgUpdaterGdal imgUpdaterEst;
754  ImgUpdaterGdal imgUpdaterUncert;
755  for(int modindex=1;modindex<nmodel;++modindex){
756  imgUpdaterEst.open(outputfw_opt[0]);
757  imgUpdaterEst.setNoData(obsnodata_opt);
758  imgUpdaterUncert.open(uncertfw_opt[0]);
759  if(verbose_opt[0]){
760  cout << "processing time " << tmodel_opt[modindex] << endl;
761  if(obsindex<relobsindex.size())
762  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
763  else
764  cout << "There is no next observation" << endl;
765  }
766 
767  //calculate regression between two subsequence model inputs
768  if(model_opt.size()==nmodel){
769  imgReaderModel1.open(model_opt[modindex-1]);
770  imgReaderModel1.setNoData(modnodata_opt);
771  imgReaderModel2.open(model_opt[modindex]);
772  imgReaderModel2.setNoData(modnodata_opt);
773  }
774  if(modelmask_opt.size()==nmodel){
775  imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
776  imgReaderModel1Mask.setNoData(msknodata_opt);
777  imgReaderModel2Mask.open(modelmask_opt[modindex]);
778  imgReaderModel2Mask.setNoData(msknodata_opt);
779  }
780 
781  pfnProgress(progress,pszMessage,pProgressArg);
782 
783  bool update=false;
784  if(obsindex<relobsindex.size()){
785  update=(relobsindex[obsindex]==modindex);
786  }
787  if(update){
788  if(observation_opt.size()==nobs){
789  if(verbose_opt[0])
790  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
791  imgReaderObs.open(observation_opt[obsindex]);
792  imgReaderObs.getGeoTransform(geotransform);
793  imgReaderObs.setNoData(obsnodata_opt);
794  }
795  if(observationmask_opt.size()==nobs){
796  imgReaderObsMask.open(observationmask_opt[obsindex]);
797  imgReaderObsMask.setNoData(msknodata_opt);
798  }
799  }
800  //prediction (also to fill cloudy pixels in measurement update mode)
801  string input;
802  input=outputfw_opt[0];
803 
804  vector< vector<double> > obsLineVector(down_opt[0]);
805  vector<double> obsLineBuffer;
806  vector<double> obsMaskLineBuffer;
807  vector<double> model1MaskLineBuffer;
808  vector<double> model2MaskLineBuffer;
809  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
810  vector<double> model1LineBuffer;
811  vector<double> model2LineBuffer;
812  vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
813  vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
814  vector<double> uncertObsLineBuffer;
815  vector< vector<double> > estLineVector(down_opt[0]);
816  vector<double> estLineBuffer;
817  vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
818  vector<double> uncertReadBuffer;
819  vector<double> estWriteBuffer(ncol);
820  vector<double> uncertWriteBuffer(ncol);
821  //test
822  vector<double> gainWriteBuffer(ncol);
823 
824  int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
825  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
826  int readModel1Band=(model_opt.size()==nmodel)? 0:modindex-1;
827  int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
828  int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex-1;
829  int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
830 
831  //initialize obsLineVector if update
832  if(update){
833  if(verbose_opt[0])
834  cout << "initialize obsLineVector" << endl;
835  assert(down_opt[0]%2);//window size must be odd
836  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
837  if(iline<0)//replicate line 0
838  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
839  else
840  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
841  }
842  }
843  //initialize estLineVector
844  if(verbose_opt[0])
845  cout << "initialize estLineVector" << endl;
846  assert(down_opt[0]%2);//window size must be odd
847 
848  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
849  if(iline<0)//replicate line 0
850  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex-1);
851  else
852  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex-1);
853  }
854  statfactory::StatFactory statobs;
855  statobs.setNoDataValues(obsnodata_opt);
856  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
857  //todo: read entire window for uncertReadBuffer...
858  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
859  imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex-1);
860  imgUpdaterUncert.image2geo(0,irow,geox,geoy);
861  if(model_opt.size()==nmodel){
862  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
863  // assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
864  imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
865  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
866  }
867  else{
868  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
869  imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
870  }
871  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
872  imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
873 
874  if(modelmask_opt.size()){
875  imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
876  if(modelmask_opt.size()==nmodel)
877  imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
878  else
879  imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
880  }
881 
882  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
883  estLineVector.erase(estLineVector.begin());
884  imgUpdaterEst.readData(estLineBuffer,maxRow,modindex-1);
885  estLineVector.push_back(estLineBuffer);
886  estLineBuffer=estLineVector[down_opt[0]/2];
887 
888  if(update){
889  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
890  obsLineVector.erase(obsLineVector.begin());
891  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
892  obsLineVector.push_back(obsLineBuffer);
893  obsLineBuffer=obsLineVector[down_opt[0]/2];
894 
895  if(observationmask_opt.size())
896  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
897  }
898 
899  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
900  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
901  imgUpdaterEst.image2geo(icol,irow,geox,geoy);
902  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
903  int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
904  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
905  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
906  estWindowBuffer.clear();
907  for(int iline=0;iline<estLineVector.size();++iline){
908  for(int isample=minCol;isample<=maxCol;++isample){
909  assert(isample<estLineVector[iline].size());
910  estWindowBuffer.push_back(estLineVector[iline][isample]);
911  }
912  }
913  if(update){
914  obsWindowBuffer.clear();
915  for(int iline=0;iline<obsLineVector.size();++iline){
916  for(int isample=minCol;isample<=maxCol;++isample){
917  assert(isample<obsLineVector[iline].size());
918  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
919  }
920  }
921  }
922  double estValue=estLineBuffer[icol];
923  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
924  bool model1IsNoData=false;
925  if(modelmask_opt.size())
926  model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
927  // lowerCol=modCol-0.5;
928  // lowerCol=static_cast<int>(lowerCol);
929  // upperCol=modCol+0.5;
930  // upperCol=static_cast<int>(upperCol);
931  lowerCol=static_cast<int>(modCol-0.5+modEps);
932  if(lowerCol<0)
933  lowerCol=0;
934  if(lowerCol>=imgReaderModel1.nrOfCol())
935  lowerCol=imgReaderModel1.nrOfCol()-1;
936  upperCol=lowerCol+1.0;
937  if(upperCol>=imgReaderModel1.nrOfCol())
938  upperCol=imgReaderModel1.nrOfCol()-1;
939  double modValue1;
940  if(!imgReaderModel1.isNoData(model1LineBuffer[lowerCol])){
941  if(!imgReaderModel1.isNoData(model1LineBuffer[upperCol])){
942  modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
943  }
944  else{
945  modValue1=model1LineBuffer[lowerCol];
946  }
947  }
948  else{
949  modValue1=model1LineBuffer[upperCol];
950  }
951  // double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
952  model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
953  if(model_opt.size()==nmodel)
954  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
955  else
956  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
957  bool model2IsNoData=false;
958  if(modelmask_opt.size())
959  model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
960  // lowerCol=modCol-0.5;
961  // lowerCol=static_cast<int>(lowerCol);
962  // upperCol=modCol+0.5;
963  // upperCol=static_cast<int>(upperCol);
964  lowerCol=static_cast<int>(modCol-0.5+modEps);
965  if(lowerCol<0)
966  lowerCol=0;
967  if(lowerCol>=imgReaderModel1.nrOfCol())
968  lowerCol=imgReaderModel1.nrOfCol()-1;
969  upperCol=lowerCol+1.0;
970  if(upperCol>=imgReaderModel1.nrOfCol())
971  upperCol=imgReaderModel1.nrOfCol()-1;
972  // double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
973  double modValue2;
974  if(!imgReaderModel1.isNoData(model2LineBuffer[lowerCol])){
975  if(!imgReaderModel1.isNoData(model2LineBuffer[upperCol])){
976  modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
977  }
978  else{
979  modValue2=model2LineBuffer[lowerCol];
980  }
981  }
982  else{
983  modValue2=model2LineBuffer[upperCol];
984  }
985  model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
986  bool obsIsNoData=false;
987  if(observationmask_opt.size())
988  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
989  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
990 
991  if(imgUpdaterEst.isNoData(estValue)){
992  //we have not found any valid data yet, better here to take the current model value if valid
993  if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
994  estWriteBuffer[icol]=obsnodata_opt[0];
995  uncertWriteBuffer[icol]=uncertNodata_opt[0];
996  //test
997  gainWriteBuffer[icol]=0;
998  }
999  else{
1000  estWriteBuffer[icol]=modValue2;
1001  uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
1002  if(obsmin_opt.size()){
1003  if(estWriteBuffer[icol]<obsmin_opt[0])
1004  estWriteBuffer[icol]=obsmin_opt[0];
1005  }
1006  if(obsmax_opt.size()){
1007  if(estWriteBuffer[icol]>obsmax_opt[0])
1008  estWriteBuffer[icol]=obsmax_opt[0];
1009  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1010  uncertWriteBuffer[icol]=obsmax_opt[0];
1011  }
1012  //test
1013  gainWriteBuffer[icol]=0;
1014  }
1015  }
1016  else{//previous estimate is valid
1017  double estMeanValue=0;
1018  double estVarValue=0;
1019  statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
1020  double nvalid=0;
1021  //time update
1022  double processNoiseVariance=processNoise_opt[0]*estVarValue;
1023  //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
1024 
1025  if(model1IsNoData||model2IsNoData){
1026  estWriteBuffer[icol]=estValue;
1027  // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
1028  //todo: check following line if makes sense
1029  uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
1030  }
1031  else{//model is good
1032  double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
1033  estWriteBuffer[icol]=estValue*modRatio;
1034  uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
1035  }
1036  if(obsmin_opt.size()){
1037  if(estWriteBuffer[icol]<obsmin_opt[0])
1038  estWriteBuffer[icol]=obsmin_opt[0];
1039  }
1040  if(obsmax_opt.size()){
1041  if(estWriteBuffer[icol]>obsmax_opt[0])
1042  estWriteBuffer[icol]=obsmax_opt[0];
1043  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1044  uncertWriteBuffer[icol]=obsmax_opt[0];
1045  }
1046  }
1047  //measurement update
1048  if(update&&!obsIsNoData){
1049  double kalmanGain=1;
1050  if(!model2IsNoData){//model is valid
1051  statfactory::StatFactory statobs;
1052  statobs.setNoDataValues(obsnodata_opt);
1053  double obsMeanValue=0;
1054  double obsVarValue=0;
1055  double difference=0;
1056  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
1057  difference=obsMeanValue-modValue2;
1058  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
1059  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
1060 
1061  if(errObs<eps_opt[0])
1062  errObs=eps_opt[0];
1063  double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
1064 
1065  if(errorCovariance+errObs>eps_opt[0])
1066  kalmanGain=errorCovariance/(errorCovariance+errObs);
1067  else
1068  kalmanGain=1;
1069  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
1070  uncertWriteBuffer[icol]*=(1-kalmanGain);
1071  if(obsmin_opt.size()){
1072  if(estWriteBuffer[icol]<obsmin_opt[0])
1073  estWriteBuffer[icol]=obsmin_opt[0];
1074  }
1075  if(obsmax_opt.size()){
1076  if(estWriteBuffer[icol]>obsmax_opt[0])
1077  estWriteBuffer[icol]=obsmax_opt[0];
1078  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1079  uncertWriteBuffer[icol]=obsmax_opt[0];
1080  }
1081  }
1082  if(kalmanGain>=1)
1083  kalmanGain=1;
1084  //test
1085  gainWriteBuffer[icol]=kalmanGain;
1086  }
1087  }
1088  }
1089 
1090  //test
1091  if(gain_opt.size())
1092  imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
1093  imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
1094  imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
1095  progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
1096  pfnProgress(progress,pszMessage,pProgressArg);
1097  }
1098  }
1099 
1100  //must close writers to ensure flush
1101  imgUpdaterEst.close();
1102  imgUpdaterUncert.close();
1103 
1104  if(update){
1105  if(observation_opt.size()==nobs)
1106  imgReaderObs.close();
1107  if(observationmask_opt.size()==nobs)
1108  imgReaderObsMask.close();
1109  ++obsindex;
1110  }
1111  if(model_opt.size()==nmodel){
1112  imgReaderModel1.close();
1113  imgReaderModel2.close();
1114  }
1115  if(modelmask_opt.size()==nmodel){
1116  imgReaderModel1Mask.close();
1117  imgReaderModel2Mask.close();
1118  }
1119  }
1120  //test
1121  if(gain_opt.size())
1122  imgWriterGain.close();
1123  }
1124  }
1125  catch(string errorString){
1126  cerr << errorString << endl;
1127  exit(1);
1128  }
1129  catch(...){
1130  cerr << "Error in forward direction " << endl;
1131  exit(2);
1132  }
1133  try{
1134  if(find(direction_opt.begin(),direction_opt.end(),"backward")!=direction_opt.end()){
1136  cout << "Running backward model" << endl;
1137  obsindex=relobsindex.size()-1;
1138  if(verbose_opt[0])
1139  cout << "Opening image " << outputbw_opt[0] << " for writing " << endl;
1140 
1141  ImgWriterGdal imgWriterEst;
1142  ImgWriterGdal imgWriterUncert;
1143  imgWriterEst.open(outputbw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1144  imgWriterEst.setProjectionProj4(projection_opt[0]);
1145  imgWriterEst.setGeoTransform(geotransform);
1146  imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
1147  imgWriterUncert.open(uncertbw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1148  imgWriterUncert.setProjectionProj4(projection_opt[0]);
1149  imgWriterUncert.setGeoTransform(geotransform);
1150 
1151  try{
1152  // //test
1153  // if(gain_opt.size()){
1154  // imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
1155  // imgWriterGain.setProjectionProj4(projection_opt[0]);
1156  // imgWriterGain.setGeoTransform(geotransform);
1157  // imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
1158  // }
1159 
1160  if(verbose_opt[0]){
1161  cout << "processing time " << tmodel_opt.back() << endl;
1162  if(obsindex<relobsindex.size()){
1163  assert(tmodel_opt.size()>relobsindex[obsindex]);
1164  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
1165  }
1166  else
1167  cout << "There is no next observation" << endl;
1168  }
1169  if(model_opt.size()==nmodel){
1170  imgReaderModel1.open(model_opt.back());
1171  imgReaderModel1.setNoData(modnodata_opt);
1172  }
1173  if(modelmask_opt.size()==nmodel){
1174  imgReaderModel1Mask.open(modelmask_opt[0]);
1175  imgReaderModel1Mask.setNoData(msknodata_opt);
1176  }
1177  }
1178  catch(string errorString){
1179  cerr << errorString << endl;
1180  }
1181  catch(...){
1182  cerr << "Error opening file " << model_opt[0] << endl;
1183  }
1184 
1185  double modEps=0.00001;
1186  double modRow=0;
1187  double modCol=0;
1188  double lowerCol=0;
1189  double upperCol=0;
1190  RESAMPLE theResample=BILINEAR;
1191 
1192  if(relobsindex.back()<nmodel-1){//initialize output_opt.back() as last model
1193  //write last model as output
1194  if(verbose_opt[0])
1195  cout << "write last model as output" << endl;
1196  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
1197  vector<double> estReadBuffer;
1198  vector<double> lineModelMask;
1199  vector<double> estWriteBuffer(ncol);
1200  vector<double> uncertWriteBuffer(ncol);
1201  // //test
1202  // vector<double> gainWriteBuffer(ncol);
1203  try{
1204  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
1205  imgWriterEst.image2geo(0,irow,geox,geoy);
1206  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1207  if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
1208  cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
1209  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1210  }
1211  // imgReaderModel1.readData(estReadBuffer,modRow,0,theResample);
1212  int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
1213  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
1214  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
1215  if(modelmask_opt.size())
1216  imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
1217  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
1218  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
1219  imgWriterEst.image2geo(icol,irow,geox,geoy);
1220  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1221  if(lineModelMask.size()>modCol){
1222  if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
1223  estWriteBuffer[icol]=obsnodata_opt[0];
1224  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1225  //test
1226  // gainWriteBuffer[icol]=obsnodata_opt[0];
1227  continue;
1228  }
1229  }
1230  // lowerCol=modCol-0.5;
1231  // lowerCol=static_cast<int>(lowerCol);
1232  // upperCol=modCol+0.5;
1233  // upperCol=static_cast<int>(upperCol);
1234  lowerCol=static_cast<int>(modCol-0.5+modEps);
1235  if(lowerCol<0)
1236  lowerCol=0;
1237  if(lowerCol>=imgReaderModel1.nrOfCol())
1238  lowerCol=imgReaderModel1.nrOfCol()-1;
1239  upperCol=lowerCol+1.0;
1240  if(upperCol>=imgReaderModel1.nrOfCol())
1241  upperCol=imgReaderModel1.nrOfCol()-1;
1242  // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
1243  double modValue;
1244  if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
1245  if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
1246  modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
1247  }
1248  else{
1249  modValue=estReadBuffer[lowerCol];
1250  }
1251  }
1252  else{
1253  modValue=estReadBuffer[upperCol];
1254  }
1255  // double modValue=estReadBuffer[modCol];
1256  if(imgReaderModel1.isNoData(modValue)){
1257  estWriteBuffer[icol]=obsnodata_opt[0];
1258  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1259  //test
1260  // gainWriteBuffer[icol]=obsnodata_opt[0];
1261  continue;
1262  }
1263  estWriteBuffer[icol]=modValue;
1264  if(obsmin_opt.size()){
1265  if(estWriteBuffer[icol]<obsmin_opt[0])
1266  estWriteBuffer[icol]=obsmin_opt[0];
1267  }
1268  if(obsmax_opt.size()){
1269  if(estWriteBuffer[icol]>obsmax_opt[0])
1270  estWriteBuffer[icol]=obsmax_opt[0];
1271  }
1272  uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
1273  //test
1274  // gainWriteBuffer[icol]=0;
1275  }
1276  }
1277  imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
1278  imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
1279  // //test
1280  // if(gain_opt.size())
1281  // imgWriterGain.writeData(gainWriteBuffer,irow,0);
1282  }
1283  }
1284  catch(string errorString){
1285  cerr << errorString << endl;
1286  }
1287  catch(...){
1288  cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
1289  }
1290  }
1291  }
1292  else{//we have a measurement at end time
1293  if(verbose_opt[0])
1294  cout << "we have a measurement at end time" << endl;
1295  if(observation_opt.size()==nobs){
1296  imgReaderObs.open(observation_opt.back());
1297  imgReaderObs.setNoData(obsnodata_opt);
1298  }
1299  if(observationmask_opt.size()==nobs){
1300  imgReaderObsMask.open(observationmask_opt.back());
1301  imgReaderObsMask.setNoData(msknodata_opt);
1302  }
1303  imgReaderObs.getGeoTransform(geotransform);
1304 
1305  vector< vector<double> > obsLineVector(down_opt[0]);
1306  vector<double> obsLineBuffer;
1307  vector<double> obsMaskLineBuffer;
1308  vector<double> modelMaskLineBuffer;
1309  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
1310  vector<double> estReadBuffer;
1311  vector<double> estWriteBuffer(ncol);
1312  vector<double> uncertWriteBuffer(ncol);
1313  vector<double> uncertObsLineBuffer;
1314  // //test
1315  // vector<double> gainWriteBuffer(ncol);
1316 
1317  if(verbose_opt[0])
1318  cout << "initialize obsLineVector" << endl;
1319  assert(down_opt[0]%2);//window size must be odd
1320  int readObsBand=(observation_opt.size()==nobs)? 0:nobs-1;
1321  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:nobs-1;
1322  int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
1323  int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:nmodel-1;
1324  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
1325  if(iline<0)//replicate line 0
1326  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
1327  else
1328  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
1329  }
1330  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
1331  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
1332  imgWriterEst.image2geo(0,irow,geox,geoy);
1333  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1334  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1335  imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
1336  if(modelmask_opt.size())
1337  imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
1338  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
1339  obsLineVector.erase(obsLineVector.begin());
1340  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
1341  obsLineVector.push_back(obsLineBuffer);
1342 
1343  if(observationmask_opt.size())
1344  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
1345 
1346  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
1347  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
1348  imgWriterEst.image2geo(icol,irow,geox,geoy);
1349  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1350  assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1351  bool modelIsNoData=false;
1352  if(modelmask_opt.size())
1353  modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
1354  // lowerCol=modCol-0.5;
1355  // lowerCol=static_cast<int>(lowerCol);
1356  // upperCol=modCol+0.5;
1357  // upperCol=static_cast<int>(upperCol);
1358  lowerCol=static_cast<int>(modCol-0.5+modEps);
1359  if(lowerCol<0)
1360  lowerCol=0;
1361  if(lowerCol>=imgReaderModel1.nrOfCol())
1362  lowerCol=imgReaderModel1.nrOfCol()-1;
1363  upperCol=lowerCol+1.0;
1364  if(upperCol>=imgReaderModel1.nrOfCol())
1365  upperCol=imgReaderModel1.nrOfCol()-1;
1366  // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
1367  double modValue;
1368  if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
1369  if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
1370  modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
1371  }
1372  else{
1373  modValue=estReadBuffer[lowerCol];
1374  }
1375  }
1376  else{
1377  modValue=estReadBuffer[upperCol];
1378  }
1379  // double modValue=estReadBuffer[modCol];
1380  double errMod=uncertModel_opt[0];//*stdDev*stdDev;
1381  modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
1382  bool obsIsNoData=false;
1383  if(observationmask_opt.size())
1384  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
1385  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
1386  if(modelIsNoData){//model is nodata: retain observation
1387  if(obsIsNoData){//both model and observation nodata
1388  estWriteBuffer[icol]=obsnodata_opt[0];
1389  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1390  //test
1391  // gainWriteBuffer[icol]=obsnodata_opt[0];
1392  }
1393  else{
1394  estWriteBuffer[icol]=obsLineBuffer[icol];
1395  if(obsmin_opt.size()){
1396  if(estWriteBuffer[icol]<obsmin_opt[0])
1397  estWriteBuffer[icol]=obsmin_opt[0];
1398  }
1399  if(obsmax_opt.size()){
1400  if(estWriteBuffer[icol]>obsmax_opt[0])
1401  estWriteBuffer[icol]=obsmax_opt[0];
1402  }
1403  uncertWriteBuffer[icol]=uncertObs_opt[0];
1404  }
1405  }
1406  else{//model is valid: calculate estimate from model
1407  estWriteBuffer[icol]=modValue;
1408  uncertWriteBuffer[icol]=errMod;//in case observation is not valid
1409  //test
1410  // gainWriteBuffer[icol]=0;
1411  }
1412  //measurement update
1413  if(!obsIsNoData){
1414  // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
1415  double kalmanGain=1;
1416  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
1417  int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
1418  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
1419  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
1420  obsWindowBuffer.clear();
1421  for(int iline=0;iline<obsLineVector.size();++iline){
1422  for(int isample=minCol;isample<=maxCol;++isample){
1423  assert(isample<obsLineVector[iline].size());
1424  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
1425  }
1426  }
1427  if(!modelIsNoData){//model is valid
1428  statfactory::StatFactory statobs;
1429  statobs.setNoDataValues(obsnodata_opt);
1430  double obsMeanValue=0;
1431  double obsVarValue=0;
1432  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
1433  double difference=0;
1434  difference=obsMeanValue-modValue;
1435  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
1436  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
1437  // double errorCovariance=errMod;
1438  double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
1439  if(errorCovariance+errObs>eps_opt[0])
1440  kalmanGain=errorCovariance/(errorCovariance+errObs);
1441  else
1442  kalmanGain=1;
1443  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
1444  if(obsmin_opt.size()){
1445  if(estWriteBuffer[icol]<obsmin_opt[0])
1446  estWriteBuffer[icol]=obsmin_opt[0];
1447  }
1448  if(obsmax_opt.size()){
1449  if(estWriteBuffer[icol]>obsmax_opt[0])
1450  estWriteBuffer[icol]=obsmax_opt[0];
1451  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1452  uncertWriteBuffer[icol]=obsmax_opt[0];
1453  }
1454  }
1455  if(kalmanGain>=1)
1456  kalmanGain=1;
1457  //test
1458  // gainWriteBuffer[icol]=kalmanGain;
1459  }
1460  }
1461  }
1462  imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
1463  imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
1464  // //test
1465  // if(gain_opt.size())
1466  // imgWriterGain.writeData(gainWriteBuffer,irow,0);
1467  }
1468  }
1469  if(observation_opt.size()==nobs)
1470  imgReaderObs.close();
1471  if(observationmask_opt.size()==nobs)
1472  imgReaderObsMask.close();
1473  --obsindex;
1474  }
1475 
1476  if(model_opt.size()==nmodel)
1477  imgReaderModel1.close();
1478  if(modelmask_opt.size()==nmodel)
1479  imgReaderModel1Mask.close();
1480  imgWriterEst.close();
1481  imgWriterUncert.close();
1482 
1483  ImgUpdaterGdal imgUpdaterEst;
1484  ImgUpdaterGdal imgUpdaterUncert;
1485  for(int modindex=nmodel-2;modindex>=0;--modindex){
1486  imgUpdaterEst.open(outputbw_opt[0]);
1487  imgUpdaterEst.setNoData(obsnodata_opt);
1488  imgUpdaterUncert.open(uncertbw_opt[0]);
1489  if(verbose_opt[0]){
1490  cout << "processing time " << tmodel_opt[modindex] << endl;
1491  if(obsindex<relobsindex.size())
1492  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
1493  else
1494  cout << "There is no next observation" << endl;
1495  }
1496 
1497  //calculate regression between two subsequence model inputs
1498  if(model_opt.size()==nmodel){
1499  imgReaderModel1.open(model_opt[modindex+1]);
1500  imgReaderModel1.setNoData(modnodata_opt);
1501  imgReaderModel2.open(model_opt[modindex]);
1502  imgReaderModel2.setNoData(modnodata_opt);
1503  }
1504  if(modelmask_opt.size()==nmodel){
1505  imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
1506  imgReaderModel1Mask.setNoData(msknodata_opt);
1507  imgReaderModel2Mask.open(modelmask_opt[modindex]);
1508  imgReaderModel2Mask.setNoData(msknodata_opt);
1509  }
1510 
1511  pfnProgress(progress,pszMessage,pProgressArg);
1512 
1513  bool update=false;
1514  if(obsindex<relobsindex.size()){
1515  update=(relobsindex[obsindex]==modindex);
1516  }
1517  if(update){
1518  if(observation_opt.size()==nobs){
1519  if(verbose_opt[0])
1520  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
1521  imgReaderObs.open(observation_opt[obsindex]);
1522  imgReaderObs.getGeoTransform(geotransform);
1523  imgReaderObs.setNoData(obsnodata_opt);
1524  }
1525  if(observationmask_opt.size()==nobs){
1526  imgReaderObsMask.open(observationmask_opt[obsindex]);
1527  imgReaderObsMask.setNoData(msknodata_opt);
1528  }
1529  }
1530  //prediction (also to fill cloudy pixels in update mode)
1531  string input;
1532  input=outputbw_opt[0];
1533 
1534  vector< vector<double> > obsLineVector(down_opt[0]);
1535  vector<double> obsLineBuffer;
1536  vector<double> obsMaskLineBuffer;
1537  vector<double> model1MaskLineBuffer;
1538  vector<double> model2MaskLineBuffer;
1539  vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
1540  vector<double> model1LineBuffer;
1541  vector<double> model2LineBuffer;
1542  vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
1543  vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
1544  vector<double> uncertObsLineBuffer;
1545  vector< vector<double> > estLineVector(down_opt[0]);
1546  vector<double> estLineBuffer;
1547  vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
1548  vector<double> uncertReadBuffer;
1549  vector<double> estWriteBuffer(ncol);
1550  vector<double> uncertWriteBuffer(ncol);
1551  //test
1552  // vector<double> gainWriteBuffer(ncol);
1553 
1554  int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
1555  int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
1556  int readModel1Band=(model_opt.size()==nmodel)? 0:modindex+1;
1557  int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
1558  int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex+1;
1559  int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
1560 
1561  //initialize obsLineVector
1562  if(update){
1563  if(verbose_opt[0])
1564  cout << "initialize obsLineVector" << endl;
1565  assert(down_opt[0]%2);//window size must be odd
1566  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
1567  if(iline<0)//replicate line 0
1568  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
1569  else
1570  imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
1571  }
1572  }
1573  //initialize estLineVector
1574  if(verbose_opt[0])
1575  cout << "initialize estLineVector" << endl;
1576  assert(down_opt[0]%2);//window size must be odd
1577 
1578  for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
1579  if(iline<0)//replicate line 0
1580  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex+1);
1581  else
1582  imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex+1);
1583  }
1584  statfactory::StatFactory statobs;
1585  statobs.setNoDataValues(obsnodata_opt);
1586 
1587  for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
1588  //todo: read entire window for uncertReadBuffer...
1589  for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
1590  imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex+1);
1591  imgUpdaterUncert.image2geo(0,irow,geox,geoy);
1592  if(model_opt.size()==nmodel){
1593  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
1594  // assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
1595  imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
1596  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1597  }
1598  else{
1599  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1600  imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
1601  }
1602 
1603  // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
1604  imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
1605  if(modelmask_opt.size()){
1606  imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
1607  if(modelmask_opt.size()==nmodel)
1608  imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
1609  else
1610  imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
1611  }
1612  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
1613  estLineVector.erase(estLineVector.begin());
1614  imgUpdaterEst.readData(estLineBuffer,maxRow,modindex+1);
1615  estLineVector.push_back(estLineBuffer);
1616  estLineBuffer=estLineVector[down_opt[0]/2];
1617 
1618  if(update){
1619  int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
1620  obsLineVector.erase(obsLineVector.begin());
1621  imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
1622  obsLineVector.push_back(obsLineBuffer);
1623  obsLineBuffer=obsLineVector[down_opt[0]/2];
1624 
1625  if(observationmask_opt.size())
1626  imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsBand);
1627  }
1628  for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
1629  for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
1630  imgUpdaterEst.image2geo(icol,irow,geox,geoy);
1631  int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
1632  int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
1633  int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
1634  int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
1635  estWindowBuffer.clear();
1636  for(int iline=0;iline<estLineVector.size();++iline){
1637  for(int isample=minCol;isample<=maxCol;++isample){
1638  assert(isample<estLineVector[iline].size());
1639  estWindowBuffer.push_back(estLineVector[iline][isample]);
1640  }
1641  }
1642  if(update){
1643  obsWindowBuffer.clear();
1644  for(int iline=0;iline<obsLineVector.size();++iline){
1645  for(int isample=minCol;isample<=maxCol;++isample){
1646  assert(isample<obsLineVector[iline].size());
1647  obsWindowBuffer.push_back(obsLineVector[iline][isample]);
1648  }
1649  }
1650  }
1651 
1652  double estValue=estLineBuffer[icol];
1653  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1654  bool model1IsNoData=false;
1655 
1656  if(modelmask_opt.size())
1657  model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
1658 
1659  // lowerCol=modCol-0.5;
1660  // lowerCol=static_cast<int>(lowerCol);
1661  // upperCol=modCol+0.5;
1662  // upperCol=static_cast<int>(upperCol);
1663  lowerCol=static_cast<int>(modCol-0.5+modEps);
1664  if(lowerCol<0)
1665  lowerCol=0;
1666  if(lowerCol>=imgReaderModel1.nrOfCol())
1667  lowerCol=imgReaderModel1.nrOfCol()-1;
1668  upperCol=lowerCol+1.0;
1669  if(upperCol>=imgReaderModel1.nrOfCol())
1670  upperCol=imgReaderModel1.nrOfCol()-1;
1671  // double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
1672  double modValue1;
1673  if(!imgReaderModel1.isNoData(model1LineBuffer[lowerCol])){
1674  if(!imgReaderModel1.isNoData(model1LineBuffer[upperCol])){
1675  modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
1676  }
1677  else{
1678  modValue1=model1LineBuffer[lowerCol];
1679  }
1680  }
1681  else{
1682  modValue1=model1LineBuffer[upperCol];
1683  }
1684  model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
1685  if(model_opt.size()==nmodel)
1686  imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
1687  else
1688  imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
1689  bool model2IsNoData=false;
1690 
1691  if(modelmask_opt.size())
1692  model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
1693  // lowerCol=modCol-0.5;
1694  // lowerCol=static_cast<int>(lowerCol);
1695  // upperCol=modCol+0.5;
1696  // upperCol=static_cast<int>(upperCol);
1697  lowerCol=static_cast<int>(modCol-0.5+modEps);
1698  if(lowerCol<0)
1699  lowerCol=0;
1700  if(lowerCol>=imgReaderModel1.nrOfCol())
1701  lowerCol=imgReaderModel1.nrOfCol()-1;
1702  upperCol=lowerCol+1.0;
1703  if(upperCol>=imgReaderModel1.nrOfCol())
1704  upperCol=imgReaderModel1.nrOfCol()-1;
1705  double modValue2;
1706  if(!imgReaderModel1.isNoData(model2LineBuffer[lowerCol])){
1707  if(!imgReaderModel1.isNoData(model2LineBuffer[upperCol])){
1708  modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
1709  }
1710  else{
1711  modValue2=model2LineBuffer[lowerCol];
1712  }
1713  }
1714  else{
1715  modValue2=model2LineBuffer[upperCol];
1716  }
1717  // double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
1718  model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
1719  bool obsIsNoData=false;
1720  if(observationmask_opt.size())
1721  obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
1722  obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
1723 
1724  if(imgUpdaterEst.isNoData(estValue)){
1725  //we have not found any valid data yet, better here to take the current model value if valid
1726  if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
1727  estWriteBuffer[icol]=obsnodata_opt[0];
1728  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1729  //test
1730  // gainWriteBuffer[icol]=0;
1731  }
1732  else{
1733  estWriteBuffer[icol]=modValue2;
1734  uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
1735  if(obsmin_opt.size()){
1736  if(estWriteBuffer[icol]<obsmin_opt[0])
1737  estWriteBuffer[icol]=obsmin_opt[0];
1738  }
1739  if(obsmax_opt.size()){
1740  if(estWriteBuffer[icol]>obsmax_opt[0])
1741  estWriteBuffer[icol]=obsmax_opt[0];
1742  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1743  uncertWriteBuffer[icol]=obsmax_opt[0];
1744  }
1745  //test
1746  // gainWriteBuffer[icol]=0;
1747  }
1748  }
1749  else{//previous estimate is valid
1750  double estMeanValue=0;
1751  double estVarValue=0;
1752  statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
1753  double nvalid=0;
1754  //time update
1755  double processNoiseVariance=processNoise_opt[0]*estVarValue;
1756  //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
1757 
1758  if(model1IsNoData||model2IsNoData){
1759  estWriteBuffer[icol]=estValue;
1760  // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
1761  //todo: check following line if makes sense
1762  uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
1763  }
1764  else{//model is good
1765  double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
1766  estWriteBuffer[icol]=estValue*modRatio;
1767  uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
1768  }
1769  if(obsmin_opt.size()){
1770  if(estWriteBuffer[icol]<obsmin_opt[0])
1771  estWriteBuffer[icol]=obsmin_opt[0];
1772  }
1773  if(obsmax_opt.size()){
1774  if(estWriteBuffer[icol]>obsmax_opt[0])
1775  estWriteBuffer[icol]=obsmax_opt[0];
1776  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1777  uncertWriteBuffer[icol]=obsmax_opt[0];
1778  }
1779  }
1780  //measurement update
1781  if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
1782  double kalmanGain=1;
1783  if(!imgReaderModel1.isNoData(modValue2)){//model is valid
1784  statfactory::StatFactory statobs;
1785  statobs.setNoDataValues(obsnodata_opt);
1786  double obsMeanValue=0;
1787  double obsVarValue=0;
1788  double difference=0;
1789  statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
1790  difference=obsMeanValue-modValue2;
1791  // errObs=uncertObs_opt[0]*sqrt(difference*difference);
1792  errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
1793 
1794  if(errObs<eps_opt[0])
1795  errObs=eps_opt[0];
1796  double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
1797 
1798  if(errorCovariance+errObs>eps_opt[0])
1799  kalmanGain=errorCovariance/(errorCovariance+errObs);
1800  else
1801  kalmanGain=1;
1802  estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
1803  uncertWriteBuffer[icol]*=(1-kalmanGain);
1804  if(obsmin_opt.size()){
1805  if(estWriteBuffer[icol]<obsmin_opt[0])
1806  estWriteBuffer[icol]=obsmin_opt[0];
1807  }
1808  if(obsmax_opt.size()){
1809  if(estWriteBuffer[icol]>obsmax_opt[0])
1810  estWriteBuffer[icol]=obsmax_opt[0];
1811  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1812  uncertWriteBuffer[icol]=obsmax_opt[0];
1813  }
1814  }
1815  if(kalmanGain>=1)
1816  kalmanGain=1;
1817  //test
1818  // gainWriteBuffer[icol]=kalmanGain;
1819  }
1820  }
1821  }
1822  // //test
1823  // if(gain_opt.size())
1824  // imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
1825  imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
1826  imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
1827  progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
1828  pfnProgress(progress,pszMessage,pProgressArg);
1829  }
1830  }
1831  //must close writers to ensure flush
1832  imgUpdaterEst.close();
1833  imgUpdaterUncert.close();
1834 
1835  if(update){
1836  if(observation_opt.size()==nobs)
1837  imgReaderObs.close();
1838  if(observationmask_opt.size()==nobs)
1839  imgReaderObsMask.close();
1840  --obsindex;
1841  }
1842  if(model_opt.size()==nmodel){
1843  imgReaderModel1.close();
1844  imgReaderModel2.close();
1845  }
1846  if(modelmask_opt.size()==nmodel){
1847  imgReaderModel1Mask.close();
1848  imgReaderModel2Mask.close();
1849  }
1850  }
1851  // //test
1852  // if(gain_opt.size())
1853  // imgWriterGain.close();
1854  }
1855  }
1856  catch(string errorString){
1857  cerr << errorString << endl;
1858  exit(1);
1859  }
1860  catch(...){
1861  cerr << "Error in backward direction " << endl;
1862  exit(2);
1863  }
1864  if(find(direction_opt.begin(),direction_opt.end(),"smooth")!=direction_opt.end()){
1866  cout << "Running smooth model" << endl;
1867  obsindex=0;
1868 
1869  ImgReaderGdal imgReaderForward(outputfw_opt[0]);
1870  ImgReaderGdal imgReaderBackward(outputbw_opt[0]);
1871  ImgReaderGdal imgReaderForwardUncert(uncertfw_opt[0]);
1872  ImgReaderGdal imgReaderBackwardUncert(uncertbw_opt[0]);
1873  imgReaderForward.setNoData(obsnodata_opt);
1874  imgReaderBackward.setNoData(obsnodata_opt);
1875 
1876  assert(imgReaderForward.nrOfBand()==nmodel);
1877  assert(imgReaderForwardUncert.nrOfBand()==nmodel);
1878  assert(imgReaderBackward.nrOfBand()==nmodel);
1879  assert(imgReaderBackwardUncert.nrOfBand()==nmodel);
1880  ImgWriterGdal imgWriterEst;
1881  imgWriterEst.setNoData(obsnodata_opt);
1882  ImgWriterGdal imgWriterUncert;
1883  imgWriterEst.open(outputfb_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1884  imgWriterEst.setProjectionProj4(projection_opt[0]);
1885  imgWriterEst.setGeoTransform(geotransform);
1886  imgWriterEst.GDALSetNoDataValue(obsnodata_opt[0]);
1887 
1888  imgWriterUncert.open(uncertfb_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
1889  imgWriterUncert.setProjectionProj4(projection_opt[0]);
1890  imgWriterUncert.setGeoTransform(geotransform);
1891  for(int modindex=0;modindex<nmodel;++modindex){
1892  if(verbose_opt[0]){
1893  cout << "processing time " << tmodel_opt[modindex] << endl;
1894  if(obsindex<relobsindex.size())
1895  cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
1896  else
1897  cout << "There is no next observation" << endl;
1898  }
1899  vector<double> estForwardBuffer;
1900  vector<double> estBackwardBuffer;
1901  vector<double> uncertObsLineBuffer;
1902  vector<double> uncertForwardBuffer;
1903  vector<double> uncertBackwardBuffer;
1904  vector<double> uncertReadBuffer;
1905  vector<double> estWriteBuffer(ncol);
1906  vector<double> uncertWriteBuffer(ncol);
1907 
1908  bool update=false;
1909  if(obsindex<relobsindex.size()){
1910  update=(relobsindex[obsindex]==modindex);
1911  }
1912 
1913  int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
1914 
1915  if(update){
1916  if(observation_opt.size()==nobs){
1917  if(verbose_opt[0])
1918  cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
1919  imgReaderObs.open(observation_opt[obsindex]);
1920  imgReaderObs.setNoData(obsnodata_opt);
1921  imgReaderObs.getGeoTransform(geotransform);
1922  }
1923  if(observationmask_opt.size()==nobs){
1924  imgReaderObsMask.open(observationmask_opt[obsindex]);
1925  imgReaderObsMask.setNoData(msknodata_opt);
1926  }
1927  }
1928 
1929  pfnProgress(progress,pszMessage,pProgressArg);
1930 
1931  for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
1932  assert(irow<imgReaderForward.nrOfRow());
1933  assert(irow<imgReaderBackward.nrOfRow());
1934  imgReaderForward.readData(estForwardBuffer,irow,modindex);
1935  imgReaderBackward.readData(estBackwardBuffer,irow,modindex);
1936  imgReaderForwardUncert.readData(uncertForwardBuffer,irow,modindex);
1937  imgReaderBackwardUncert.readData(uncertBackwardBuffer,irow,modindex);
1938 
1939  if(update){
1940  if(observation_opt.size()==nobs)
1941  imgReaderObs.readData(estWriteBuffer,irow,readObsBand);
1942  if(observationmask_opt.size())
1943  imgReaderObsMask.readData(uncertObsLineBuffer,irow,readObsBand);
1944  }
1945 
1946  for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
1947  imgWriterEst.image2geo(icol,irow,geox,geoy);
1948  double A=estForwardBuffer[icol];
1949  double B=estBackwardBuffer[icol];
1950  double C=uncertForwardBuffer[icol];
1951  double D=uncertBackwardBuffer[icol];
1952  double uncertObs=uncertObs_opt[0];
1953 
1954  double noemer=(C+D);
1955  //todo: consistently check for division by zero...
1956  if(imgReaderForward.isNoData(A)&&imgReaderBackward.isNoData(B)){
1957  estWriteBuffer[icol]=obsnodata_opt[0];
1958  uncertWriteBuffer[icol]=uncertNodata_opt[0];
1959  }
1960  else if(imgReaderForward.isNoData(A)){
1961  estWriteBuffer[icol]=B;
1962  uncertWriteBuffer[icol]=uncertBackwardBuffer[icol];
1963  }
1964  else if(imgReaderForward.isNoData(B)){
1965  estWriteBuffer[icol]=A;
1966  uncertWriteBuffer[icol]=uncertForwardBuffer[icol];
1967  }
1968  else{
1969  if(noemer<eps_opt[0]){//simple average if both uncertainties are ~>0
1970  estWriteBuffer[icol]=0.5*(A+B);
1971  uncertWriteBuffer[icol]=eps_opt[0];
1972  }
1973  else{
1974  estWriteBuffer[icol]=(A*D+B*C)/noemer;
1975  uncertWriteBuffer[icol]=C*D/noemer;
1976  if(obsmin_opt.size()){
1977  if(estWriteBuffer[icol]<obsmin_opt[0])
1978  estWriteBuffer[icol]=obsmin_opt[0];
1979  }
1980  if(obsmax_opt.size()){
1981  if(estWriteBuffer[icol]>obsmax_opt[0])
1982  estWriteBuffer[icol]=obsmax_opt[0];
1983  if(uncertWriteBuffer[icol]>obsmax_opt[0])
1984  uncertWriteBuffer[icol]=obsmax_opt[0];
1985  }
1986  }
1987  }
1988  }
1989  imgWriterEst.writeData(estWriteBuffer,irow,modindex);
1990  imgWriterUncert.writeData(uncertWriteBuffer,irow,modindex);
1991  progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
1992  pfnProgress(progress,pszMessage,pProgressArg);
1993  }
1994  if(update){
1995  if(observation_opt.size()==nobs)
1996  imgReaderObs.close();
1997  ++obsindex;
1998  }
1999  }
2000  imgReaderForward.close();
2001  imgReaderBackward.close();
2002  imgWriterEst.close();
2003  imgWriterUncert.close();
2004  }
2005  if(observation_opt.size()<nobs)
2006  imgReaderObs.close();
2007  if(model_opt.size()<nmodel)
2008  imgReaderModel1.close();
2009 }
void close(void)
Close the image.
STL namespace.
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row) ...
bool isNoData(double value) const
Check if value is nodata in this dataset.
CPLErr setProjectionProj4(const std::string &projection)
Set the projection for this dataset from user input (supports epsg:<number> format) ...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
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
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) ...
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
double getDeltaX(void) const
Get the pixel cell spacing in x.
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...
std::string getFileName() const
Get the filename of this dataset.
Definition: ImgRasterGdal.h:96
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
void close(void)
Close the image.
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
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 nrOfBand(void) const
Get the number of bands of this dataset.
std::string getInterleave() const
Get the band coding (interleave)