| 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
259 statfactory::StatFactory stat;
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}
ImgRasterGdal::setNoData
int setNoData(const std::vector< double > nodata)
Set the no data values of this dataset using a standard template library (stl) vector as input.
Definition: ImgRasterGdal.h:140
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::geo2image
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
Definition: ImgRasterGdal.cc:387
ImgRasterGdal::getProjection
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
Definition: ImgRasterGdal.cc:40
ImgRasterGdal::isNoData
bool isNoData(double value) const
Check if value is nodata in this dataset.
Definition: ImgRasterGdal.h:136
ImgRasterGdal::getGeoTransform
std::string getGeoTransform() const
Get the geotransform data for this dataset as a string.
Definition: ImgRasterGdal.cc:200
ImgRasterGdal::setGeoTransform
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
Definition: ImgRasterGdal.cc:148
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::setProjectionProj4
CPLErr setProjectionProj4(const std::string &projection)
Set the projection for this dataset from user input (supports epsg:<number> format)
Definition: ImgRasterGdal.cc:63
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgRasterGdal::getDeltaX
double getDeltaX(void) const
Get the pixel cell spacing in x.
Definition: ImgRasterGdal.h:152
ImgRasterGdal::getInterleave
std::string getInterleave() const
Get the band coding (interleave)
Definition: ImgRasterGdal.cc:293
ImgRasterGdal::image2geo
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
Definition: ImgRasterGdal.cc:416
ImgRasterGdal::GDALSetNoDataValue
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...
Definition: ImgRasterGdal.h:142
ImgRasterGdal::getDataType
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
Definition: ImgRasterGdal.cc:92
ImgRasterGdal::getFileName
std::string getFileName() const
Get the filename of this dataset.
Definition: ImgRasterGdal.h:96
ImgReaderGdal
Definition: ImgReaderGdal.h:37
ImgReaderGdal::readData
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
ImgReaderGdal::close
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
Definition: ImgReaderGdal.cc:42
ImgReaderGdal::open
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
Definition: ImgReaderGdal.cc:36
ImgUpdaterGdal
Definition: ImgGdal.h:28
ImgUpdaterGdal::close
void close(void)
Close the image.
Definition: ImgUpdaterGdal.cc:51
ImgWriterGdal
Definition: ImgWriterGdal.h:38
ImgWriterGdal::close
void close(void)
Close the image.
Definition: ImgWriterGdal.cc:74
ImgWriterGdal::open
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...
Definition: ImgWriterGdal.cc:43
ImgWriterGdal::writeData
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
Optionpk
Definition: Optionpk.h:107
statfactory::StatFactory
Definition: StatFactory.h:43