[pktools] 02/09: New upstream version 2.6.7.2+ds
Bas Couwenberg
sebastic at debian.org
Wed Feb 7 09:03:10 UTC 2018
This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository pktools.
commit 2116a874d4a71dfb22bbcbc4a1163ae4b92d867d
Author: Bas Couwenberg <sebastic at xs4all.nl>
Date: Wed Feb 7 07:59:16 2018 +0100
New upstream version 2.6.7.2+ds
---
CMakeFiles/cmake.check_cache | 1 +
CMakeLists.txt | 4 +-
src/algorithms/Filter.cc | 213 +--
src/apps/pkcomposite.cc | 655 ++++-----
src/apps/pkcrop.cc | 538 ++++---
src/apps/pkkalman.cc | 2873 ++++++++++++++++++++------------------
src/apps/pkstat.cc | 858 ++++++------
src/apps/pkstatprofile.cc | 4 +-
src/imageclasses/ImgReaderGdal.h | 38 +-
9 files changed, 2676 insertions(+), 2508 deletions(-)
diff --git a/CMakeFiles/cmake.check_cache b/CMakeFiles/cmake.check_cache
new file mode 100644
index 0000000..3dccd73
--- /dev/null
+++ b/CMakeFiles/cmake.check_cache
@@ -0,0 +1 @@
+# This file is generated by cmake for dependency checking of the CMakeCache.txt file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 79cc307..6d730fe 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,7 +35,7 @@ set (PROJECT_SOURCE_DIR src)
# The version number.
set (PKTOOLS_VERSION_MAJOR 2)
set (PKTOOLS_VERSION_MINOR 6)
-set (PKTOOLS_VERSION_PATCH 7.1)
+set (PKTOOLS_VERSION_PATCH 7.2)
set (PKTOOLS_VERSION "${PKTOOLS_VERSION_MAJOR}.${PKTOOLS_VERSION_MINOR}.${PKTOOLS_VERSION_PATCH}")
set (PKTOOLS_PACKAGE_VERSION "${PKTOOLS_VERSION_MAJOR}.${PKTOOLS_VERSION_MINOR}.${PKTOOLS_VERSION_PATCH}")
set (PKTOOLS_PACKAGE_STRING "PKTOOLS ${PKTOOLS_VERSION_MAJOR}.${PKTOOLS_VERSION_MINOR}.${PKTOOLS_VERSION_PATCH}")
@@ -398,11 +398,13 @@ list(APPEND CPACK_SOURCE_IGNORE_FILES ".gz")
list(APPEND CPACK_SOURCE_IGNORE_FILES ".bz2")
list(APPEND CPACK_SOURCE_IGNORE_FILES ".zip")
list(APPEND CPACK_SOURCE_IGNORE_FILES ".svn")
+list(APPEND CPACK_SOURCE_IGNORE_FILES ".git")
list(APPEND CPACK_SOURCE_IGNORE_FILES "README")
list(APPEND CPACK_SOURCE_IGNORE_FILES "HOWTORELEASE.txt")
list(APPEND CPACK_SOURCE_IGNORE_FILES "CMakeCache.txt")
list(APPEND CPACK_SOURCE_IGNORE_FILES "CPackConfig.cmake")
list(APPEND CPACK_SOURCE_IGNORE_FILES "schemas")
+list(APPEND CPACK_SOURCE_IGNORE_FILES "/build/;~$;${CPACK_SOURCE_IGNORE_FILES}")
include(CPack)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index c7d954c..8bc8f03 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -72,7 +72,7 @@ void filter::Filter::dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, con
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
@@ -82,7 +82,7 @@ void filter::Filter::dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, con
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -103,7 +103,7 @@ void filter::Filter::dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, con
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
@@ -113,7 +113,7 @@ void filter::Filter::dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, con
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -134,7 +134,7 @@ void filter::Filter::dwtCut(ImgReaderGdal& input, ImgWriterGdal& output, const s
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
@@ -144,7 +144,7 @@ void filter::Filter::dwtCut(ImgReaderGdal& input, ImgWriterGdal& output, const s
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -165,22 +165,22 @@ void filter::Filter::dwtCutFrom(ImgReaderGdal& input, ImgWriterGdal& output, con
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
dwtForward(pixelInput,wavelet_type,family);
for(int iband=0;iband<input.nrOfBand();++iband){
- if(iband>=band)
- pixelInput[iband]=0;
+ if(iband>=band)
+ pixelInput[iband]=0;
}
dwtInverse(pixelInput,wavelet_type,family);
for(int iband=0;iband<input.nrOfBand();++iband)
- lineOutput[iband][x]=pixelInput[iband];
+ lineOutput[iband][x]=pixelInput[iband];
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -197,7 +197,7 @@ void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wa
//make sure data size if power of 2
while(data.size()&(data.size()-1))
data.push_back(data.back());
-
+
int nsize=data.size();
gsl_wavelet *w;
gsl_wavelet_workspace *work;
@@ -269,7 +269,7 @@ void filter::Filter::morphology(ImgReaderGdal& input, ImgWriterGdal& output, con
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
@@ -281,7 +281,7 @@ void filter::Filter::morphology(ImgReaderGdal& input, ImgWriterGdal& output, con
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -303,7 +303,7 @@ void filter::Filter::smoothNoData(ImgReaderGdal& input, const std::string& inter
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
@@ -314,7 +314,7 @@ void filter::Filter::smoothNoData(ImgReaderGdal& input, const std::string& inter
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -345,7 +345,7 @@ void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output)
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
vector<double> pixelOutput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
@@ -356,7 +356,7 @@ void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output)
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -381,44 +381,44 @@ void filter::Filter::stat(ImgReaderGdal& input, ImgWriterGdal& output, const std
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
pixelInput=lineInput.selectCol(x);
switch(getFilterType(method)){
case(filter::median):
- lineOutput[x]=stat.median(pixelInput);
- break;
+ lineOutput[x]=stat.median(pixelInput);
+ break;
case(filter::min):
- lineOutput[x]=stat.mymin(pixelInput);
- break;
+ lineOutput[x]=stat.mymin(pixelInput);
+ break;
case(filter::max):
- lineOutput[x]=stat.mymax(pixelInput);
- break;
+ lineOutput[x]=stat.mymax(pixelInput);
+ break;
case(filter::sum):
- lineOutput[x]=stat.sum(pixelInput);
- break;
+ lineOutput[x]=stat.sum(pixelInput);
+ break;
case(filter::var):
- lineOutput[x]=stat.var(pixelInput);
- break;
+ lineOutput[x]=stat.var(pixelInput);
+ break;
case(filter::stdev):
- lineOutput[x]=sqrt(stat.var(pixelInput));
- break;
+ lineOutput[x]=sqrt(stat.var(pixelInput));
+ break;
case(filter::mean):
- lineOutput[x]=stat.mean(pixelInput);
- break;
+ lineOutput[x]=stat.mean(pixelInput);
+ break;
case(filter::percentile):
- assert(m_threshold.size());
- lineOutput[x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),m_threshold[0]);
- break;
+ assert(m_threshold.size());
+ lineOutput[x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),m_threshold[0]);
+ break;
default:
- std::string errorString="method not supported";
- throw(errorString);
- break;
+ std::string errorString="method not supported";
+ throw(errorString);
+ break;
}
}
try{
- output.writeData(lineOutput,GDT_Float64,y);
+ output.writeData(lineOutput,y);
}
catch(string errorstring){
cerr << errorstring << "in line " << y << endl;
@@ -443,59 +443,64 @@ void filter::Filter::stats(ImgReaderGdal& input, ImgWriterGdal& output, const ve
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
for(int x=0;x<input.nrOfCol();++x){
- pixelInput=lineInput.selectCol(x);
- int ithreshold=0;//threshold to use for percentiles
- for(int imethod=0;imethod<methods.size();++imethod){
- switch(getFilterType(methods[imethod])){
- case(filter::nvalid):
- lineOutput[imethod][x]=stat.nvalid(pixelInput);
- break;
- case(filter::median):
- lineOutput[imethod][x]=stat.median(pixelInput);
- break;
- case(filter::min):
- lineOutput[imethod][x]=stat.mymin(pixelInput);
- break;
- case(filter::max):
- lineOutput[imethod][x]=stat.mymax(pixelInput);
- break;
- case(filter::sum):
- lineOutput[imethod][x]=stat.sum(pixelInput);
- break;
- case(filter::var):
- lineOutput[imethod][x]=stat.var(pixelInput);
- break;
- case(filter::stdev):
- lineOutput[imethod][x]=sqrt(stat.var(pixelInput));
- break;
- case(filter::mean):
- lineOutput[imethod][x]=stat.mean(pixelInput);
- break;
- case(filter::percentile):{
- assert(m_threshold.size());
- double threshold=(ithreshold<m_threshold.size())? m_threshold[ithreshold] : m_threshold[0];
- lineOutput[imethod][x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),threshold);
- ++ithreshold;
- break;
- }
- default:
- std::string errorString="method not supported";
- throw(errorString);
- break;
- }
+ try{
+ pixelInput=lineInput.selectCol(x);
}
- }
- for(int imethod=0;imethod<methods.size();++imethod){
+ catch(...){
+ std::cout << "error is caught" << std::endl;
+ exit(1);
+ }
+ int ithreshold=0;//threshold to use for percentiles
try{
- output.writeData(lineOutput[imethod],GDT_Float64,y,imethod);
+ for(int imethod=0;imethod<methods.size();++imethod){
+ switch(getFilterType(methods[imethod])){
+ case(filter::nvalid):
+ lineOutput[imethod][x]=stat.nvalid(pixelInput);
+ break;
+ case(filter::median):
+ lineOutput[imethod][x]=stat.median(pixelInput);
+ break;
+ case(filter::min):
+ lineOutput[imethod][x]=stat.mymin(pixelInput);
+ break;
+ case(filter::max):
+ lineOutput[imethod][x]=stat.mymax(pixelInput);
+ break;
+ case(filter::sum):
+ lineOutput[imethod][x]=stat.sum(pixelInput);
+ break;
+ case(filter::var):
+ lineOutput[imethod][x]=stat.var(pixelInput);
+ break;
+ case(filter::stdev):
+ lineOutput[imethod][x]=sqrt(stat.var(pixelInput));
+ break;
+ case(filter::mean):
+ lineOutput[imethod][x]=stat.mean(pixelInput);
+ break;
+ case(filter::percentile):{
+ assert(m_threshold.size());
+ double threshold=(ithreshold<m_threshold.size())? m_threshold[ithreshold] : m_threshold[0];
+ lineOutput[imethod][x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),threshold);
+ ++ithreshold;
+ break;
+ }
+ default:
+ std::string errorString="method not supported";
+ throw(errorString);
+ break;
+ }
+ }
}
- catch(string errorstring){
- cerr << errorstring << "in line " << y << endl;
+ catch(...){
+ cerr << "An Error in line " << y << endl;
}
}
+ for(int imethod=0;imethod<methods.size();++imethod)
+ output.writeData(lineOutput[imethod],y,imethod);
progress=(1.0+y)/output.nrOfRow();
pfnProgress(progress,pszMessage,pProgressArg);
}
@@ -512,7 +517,7 @@ void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output, const s
pfnProgress(progress,pszMessage,pProgressArg);
for(int y=0;y<input.nrOfRow();++y){
for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
+ input.readData(lineInput[iband],y,iband);
vector<double> pixelInput(input.nrOfBand());
vector<double> pixelOutput;
for(int x=0;x<input.nrOfCol();++x){
@@ -520,13 +525,13 @@ void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output, const s
filter(pixelInput,pixelOutput,method,dim);
for(int iband=0;iband<pixelOutput.size();++iband){
lineOutput[iband][x]=pixelOutput[iband];
- // if(pixelInput[iband]!=0)
- // assert(pixelOutput[iband]!=0);
+ // if(pixelInput[iband]!=0)
+ // assert(pixelOutput[iband]!=0);
}
}
for(int iband=0;iband<input.nrOfBand();++iband){
try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ output.writeData(lineOutput[iband],y,iband);
}
catch(string errorstring){
cerr << errorstring << "in band " << iband << ", line " << y << endl;
@@ -571,7 +576,7 @@ void filter::Filter::getSavGolayCoefficients(vector<double> &tapz, int np, int n
// for(kk = 0; kk < np; ++kk)
// c[kk] = 0.0;
for(k = -nl; k <= nr; ++k) {
- // for(k = -nl; k < nr; ++k) {
+ // for(k = -nl; k < nr; ++k) {
sum = b[0];
fac = 1.0;
for(mm = 1; mm <= m; ++mm)
@@ -618,26 +623,26 @@ void filter::Filter::ludcmp(vector<double> &a, vector<int> &indx, double &d) {
for(i = 0; i < j; ++i) {
sum = a[i * n + j];
for(k = 0; k < i; ++k)
- sum -= a[i * n + k] * a[k * n + j];
+ sum -= a[i * n + k] * a[k * n + j];
a[i * n + j] = sum;
}
big = 0.0;
for(i = j; i < n; ++i) {
sum = a[i * n + j];
for(k = 0; k < j; ++k)
- sum -= a[i * n + k] * a[k * n + j];
+ sum -= a[i * n + k] * a[k * n + j];
a[i * n + j] = sum;
if((dum = vv[i] * fabs(sum)) >= big) {
- big = dum;
- imax = i;
+ big = dum;
+ imax = i;
}
}
if(j != imax) {
for(k = 0; k < n; ++k) {
- dum = a[imax * n + k];
- a[imax * n + k] = a[j * n + k];
- a[j * n + k] = dum;
+ dum = a[imax * n + k];
+ a[imax * n + k] = a[j * n + k];
+ a[j * n + k] = dum;
}
d = -d;
vv[imax] = vv[j];
@@ -647,7 +652,7 @@ void filter::Filter::ludcmp(vector<double> &a, vector<int> &indx, double &d) {
if(j != n - 1) {
dum = 1. / a[j * n + j];
for(i = j + 1; i < n; ++i)
- a[i * n + j] *= dum;
+ a[i * n + j] *= dum;
}
}
}
@@ -663,7 +668,7 @@ void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double>
b[ip] = b[i];
if(ii != 0)
for(j = ii - 1; j < i; ++j)
- sum -= a[i * n + j] * b[j];
+ sum -= a[i * n + j] * b[j];
else if(sum != 0.0)
ii = i + 1;
b[i] = sum;
@@ -677,9 +682,9 @@ void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double>
}
double filter::Filter::getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta, bool verbose)
-{
+{
assert(srf.size()==2);//[0]: wavelength, [1]: response function
- int nband=srf[0].size();
+ int nband=srf[0].size();
double start=floor(wavelengthIn[0]);
double end=ceil(wavelengthIn.back());
if(verbose)
@@ -699,7 +704,7 @@ double filter::Filter::getCentreWavelength(const std::vector<double> &wavelength
if(verbose)
std::cout << "norm of srf: " << norm << std::endl << std::flush;
gsl_spline_free(spline);
- gsl_interp_accel_free(acc);
+ gsl_interp_accel_free(acc);
std::vector<double> wavelength_fine;
for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
wavelength_fine.push_back(win);
@@ -724,7 +729,7 @@ double filter::Filter::getCentreWavelength(const std::vector<double> &wavelength
stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
-
+
gsl_spline_free(splineOut);
gsl_interp_accel_free(accOut);
diff --git a/src/apps/pkcomposite.cc b/src/apps/pkcomposite.cc
index 4dcc156..37f8030 100644
--- a/src/apps/pkcomposite.cc
+++ b/src/apps/pkcomposite.cc
@@ -21,6 +21,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <vector>
#include <iostream>
#include <string>
+#include <math.h>
#include "imageclasses/ImgReaderGdal.h"
#include "imageclasses/ImgWriterGdal.h"
#include "imageclasses/ImgReaderOgr.h"
@@ -31,108 +32,108 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
/******************************************************************************/
/*! \page pkcomposite pkcomposite
- program to mosaic and composite geo-referenced images
-## SYNOPSIS
+ program to mosaic and composite geo-referenced images
+ ## SYNOPSIS
-<code>
+ <code>
Usage: pkcomposite -i input [-i input]* -o output
-</code>
+ </code>
-<code>
+ <code>
Options: [-b band]* [-dx xres] [-dy yres] [-e vector] [-ulx ULX -uly ULY -lrx LRX -lry LRY] [-cr rule] [-cb band] [-srcnodata value] [-bndnodata band] [-min value] [-max value] [-dstnodata value] [-r resampling_method] [-ot {Byte / Int16 / UInt16 / UInt32 / Int32 / Float32 / Float64 / CInt16 / CInt32 / CFloat32 / CFloat64}] [-of format] [-co NAME=VALUE]* [-a_srs epsg:number]
Advanced options:
- [-file] [-w weight]* [-c class]* [-ct colortable] [-d description] [-align]
-</code>
+ [-file] [-w weight]* [-c class]* [-ct colortable] [-d description] [-align]
+ </code>
-\section pkcomposite_description Description
+ \section pkcomposite_description Description
-The utility pkcomposite can be used to \em mosaic and \em composite multiple (georeferenced) raster datasets. A mosaic can merge images with different geographical extents into a single larger image. Compositing resolves the overlapping pixels according to some rule (e.g, the median of all overlapping pixels). This utility is complementary to GDAL, which currently does not support a composite step. Input datasets can have different bounding boxes and spatial resolutionsresolution.
+ The utility pkcomposite can be used to \em mosaic and \em composite multiple (georeferenced) raster datasets. A mosaic can merge images with different geographical extents into a single larger image. Compositing resolves the overlapping pixels according to some rule (e.g, the median of all overlapping pixels). This utility is complementary to GDAL, which currently does not support a composite step. Input datasets can have different bounding boxes and spatial resolutionsresolution.
-\anchor pkcomposite_rules
-composite rule | composite output
-------------- | -------------
-overwrite | Overwrite overlapping pixels: the latter input image on the command line overrules the previous image
-maxndvi | Create a maximum NDVI (normalized difference vegetation index) composite from multi-band input images. Use option -cb to set the indexes of the red and near infrared bands respectively (e.g., -cb 0 -cb 1)
-maxband | Select the pixel with a maximum value in the band specified by option -cb
-minband | Select the pixel with a minimum value in the band specified by option -cb
-mean | Calculate the mean (average) of overlapping pixels
-stdev | Calculate the standard deviation of overlapping pixels
-median | Calculate the median of overlapping pixels
-mode | Select the mode of overlapping pixels (maximum voting): use for Byte images only
-sum | Calculate the arithmetic sum of overlapping pixels
-maxallbands | For each individual band, assign the maximum value found in all overlapping pixels. Unlike maxband, output band values cannot be attributed to a single (date) pixel in the input time series
-minallbands | For each individual band, assign the minimum value found in all overlapping pixels. Unlike minband, output band values cannot be attributed to a single (date) pixel in the input time series
+ \anchor pkcomposite_rules
+ composite rule | composite output
+ ------------- | -------------
+ overwrite | Overwrite overlapping pixels: the latter input image on the command line overrules the previous image
+ maxndvi | Create a maximum NDVI (normalized difference vegetation index) composite from multi-band input images. Use option -cb to set the indexes of the red and near infrared bands respectively (e.g., -cb 0 -cb 1)
+ maxband | Select the pixel with a maximum value in the band specified by option -cb
+ minband | Select the pixel with a minimum value in the band specified by option -cb
+ mean | Calculate the mean (average) of overlapping pixels
+ stdev | Calculate the standard deviation of overlapping pixels
+ median | Calculate the median of overlapping pixels
+ mode | Select the mode of overlapping pixels (maximum voting): use for Byte images only
+ sum | Calculate the arithmetic sum of overlapping pixels
+ maxallbands | For each individual band, assign the maximum value found in all overlapping pixels. Unlike maxband, output band values cannot be attributed to a single (date) pixel in the input time series
+ minallbands | For each individual band, assign the minimum value found in all overlapping pixels. Unlike minband, output band values cannot be attributed to a single (date) pixel in the input time series
-Example: Calculate the maximum NDVI composite of two multispectral input images (e.g., red is band 0 and near infrared is band 1)
+ Example: Calculate the maximum NDVI composite of two multispectral input images (e.g., red is band 0 and near infrared is band 1)
-\code
-pkcomposite -i input1.tif -i input2.tif -o output.tif -cr maxndvi -cb 0 -cb 1
-\endcode
+ \code
+ pkcomposite -i input1.tif -i input2.tif -o output.tif -cr maxndvi -cb 0 -cb 1
+ \endcode
-Example: Calculate the minimum nadir composite of two input images, where the forth band (b=3) contains the view zenith angle
+ Example: Calculate the minimum nadir composite of two input images, where the forth band (b=3) contains the view zenith angle
-\code
-pkcomposite -i input1.tif -i input2.tif -o minzenith.tif -cr minband -cb 3
-\endcode
+ \code
+ pkcomposite -i input1.tif -i input2.tif -o minzenith.tif -cr minband -cb 3
+ \endcode
-Example: Calculate the minimum of two input images in all bands
+ Example: Calculate the minimum of two input images in all bands
-\code
-pkcomposite -i input1.tif -i input2.tif -o minimum.tif -cr minallbands
-\endcode
+ \code
+ pkcomposite -i input1.tif -i input2.tif -o minimum.tif -cr minallbands
+ \endcode
-\section pkcomposite_options Options
- - use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
- - short option `-h` shows basic options only, long option `--help` shows all options
+ \section pkcomposite_options Options
+ - use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
+ - short option `-h` shows basic options only, long option `--help` shows all options
-|short|long|type|default|description|
-|-----|----|----|-------|-----------|
- | i | input | std::string | |Input image file(s). If input contains multiple images, a multi-band output is created |
- | o | output | std::string | |Output image file |
- | b | band | int | |band index(es) to crop (leave empty if all bands must be retained) |
- | dx | dx | double | |Output resolution in x (in meter) (empty: keep original resolution) |
- | dy | dy | double | |Output resolution in y (in meter) (empty: keep original resolution) |
- | e | extent | std::string | |get boundary from extent from polygons in vector file |
- | cut | crop_to_cutline | bool | false |Crop the extent of the target dataset to the extent of the cutline |
- | eo | eo | std::string | |special extent options controlling rasterization: ATTRIBUTE or CHUNKYSIZE or ALL_TOUCHED or BURN_VALUE_FROM or MERGE_ALG, e.g., -eo ATTRIBUTE=fieldname |
- | m | mask | std::string | |Use the first band of the specified file as a validity mask (0 is nodata) |
- | msknodata | msknodata | float | 0 |Mask value not to consider for composite
- | mskband | mskband | short | 0 |Mask band to read (0 indexed) |
- | ulx | ulx | double | 0 |Upper left x value bounding box |
- | uly | uly | double | 0 |Upper left y value bounding box |
- | lrx | lrx | double | 0 |Lower right x value bounding box |
- | lry | lry | double | 0 |Lower right y value bounding box |
- | cr | crule | std::string | overwrite |Composite rule (overwrite, maxndvi, maxband, minband, mean, mode (only for byte images), median, sum, maxallbands, minallbands, stdev |
- | cb | cband | int | 0 |band index used for the composite rule (e.g., for ndvi, use --cband=0 --cband=1 with 0 and 1 indices for red and nir band respectively |
- | srcnodata | srcnodata | double | |invalid value(s) for input raster dataset |
- | bndnodata | bndnodata | int | 0 |Band(s) in input image to check if pixel is valid (used for srcnodata, min and max options) |
- | min | min | double | |flag values smaller or equal to this value as invalid. |
- | max | max | double | |flag values larger or equal to this value as invalid. |
- | dstnodata | dstnodata | double | 0 |nodata value to put in output raster dataset if not valid or out of bounds. |
- | r | resampling-method | std::string | near |Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation). |
- | ot | otype | std::string | |Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image |
- | of | oformat | std::string | GTiff |Output image format (see also gdal_translate).|
- | co | co | std::string | |Creation option for output file. Multiple options can be specified. |
- | a_srs | a_srs | std::string | |Override the spatial reference for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid |
- | file | file | short | 0 |write number of observations (1) or sequence nr of selected file (2) for each pixels as additional layer in composite |
- | w | weight | short | 1 |Weights (type: short) for the composite, use one weight for each input file in same order as input files are provided). Use value 1 for equal weights. |
- | c | class | short | 0 |classes for multi-band output image: each band represents the number of observations for one specific class. Use value 0 for no multi-band output image. |
- | ct | ct | std::string | |color table file with 5 columns: id R G B ALFA (0: transparent, 255: solid) |
- | align | align | bool | |Align output bounding box to first input image |
- | scale | scale | double | |output=scale*input+offset |
- | off | offset | double | |output=scale*input+offset |
- | d | description | std::string | |Set image description |
+ |short|long|type|default|description|
+ |-----|----|----|-------|-----------|
+ | i | input | std::string | |Input image file(s). If input contains multiple images, a multi-band output is created |
+ | o | output | std::string | |Output image file |
+ | b | band | int | |band index(es) to crop (leave empty if all bands must be retained) |
+ | dx | dx | double | |Output resolution in x (in meter) (empty: keep original resolution) |
+ | dy | dy | double | |Output resolution in y (in meter) (empty: keep original resolution) |
+ | e | extent | std::string | |get boundary from extent from polygons in vector file |
+ | cut | crop_to_cutline | bool | false |Crop the extent of the target dataset to the extent of the cutline |
+ | eo | eo | std::string | |special extent options controlling rasterization: ATTRIBUTE or CHUNKYSIZE or ALL_TOUCHED or BURN_VALUE_FROM or MERGE_ALG, e.g., -eo ATTRIBUTE=fieldname |
+ | m | mask | std::string | |Use the first band of the specified file as a validity mask (0 is nodata) |
+ | msknodata | msknodata | float | 0 |Mask value not to consider for composite
+ | mskband | mskband | short | 0 |Mask band to read (0 indexed) |
+ | ulx | ulx | double | 0 |Upper left x value bounding box |
+ | uly | uly | double | 0 |Upper left y value bounding box |
+ | lrx | lrx | double | 0 |Lower right x value bounding box |
+ | lry | lry | double | 0 |Lower right y value bounding box |
+ | cr | crule | std::string | overwrite |Composite rule (overwrite, maxndvi, maxband, minband, mean, mode (only for byte images), median, sum, maxallbands, minallbands, stdev |
+ | cb | cband | int | 0 |band index used for the composite rule (e.g., for ndvi, use --cband=0 --cband=1 with 0 and 1 indices for red and nir band respectively |
+ | srcnodata | srcnodata | double | |invalid value(s) for input raster dataset |
+ | bndnodata | bndnodata | int | 0 |Band(s) in input image to check if pixel is valid (used for srcnodata, min and max options) |
+ | min | min | double | |flag values smaller or equal to this value as invalid. |
+ | max | max | double | |flag values larger or equal to this value as invalid. |
+ | dstnodata | dstnodata | double | 0 |nodata value to put in output raster dataset if not valid or out of bounds. |
+ | r | resampling-method | std::string | near |Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation). |
+ | ot | otype | std::string | |Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image |
+ | of | oformat | std::string | GTiff |Output image format (see also gdal_translate).|
+ | co | co | std::string | |Creation option for output file. Multiple options can be specified. |
+ | a_srs | a_srs | std::string | |Override the spatial reference for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid |
+ | file | file | short | 0 |write number of observations (1) or sequence nr of selected file (2) for each pixels as additional layer in composite |
+ | w | weight | short | 1 |Weights (type: short) for the composite, use one weight for each input file in same order as input files are provided). Use value 1 for equal weights. |
+ | c | class | short | 0 |classes for multi-band output image: each band represents the number of observations for one specific class. Use value 0 for no multi-band output image. |
+ | ct | ct | std::string | |color table file with 5 columns: id R G B ALFA (0: transparent, 255: solid) |
+ | align | align | bool | |Align output bounding box to first input image |
+ | scale | scale | double | |output=scale*input+offset |
+ | off | offset | double | |output=scale*input+offset |
+ | d | description | std::string | |Set image description |
-Examples
-========
-Some examples how to use pkcomposite can be found \ref examples_pkcomposite "here"
+ Examples
+ ========
+ Some examples how to use pkcomposite can be found \ref examples_pkcomposite "here"
-FAQ
-===
-Frequently asked questions on pkcomposite can be found \ref faq_pkcomposite "here
+ FAQ
+ ===
+ Frequently asked questions on pkcomposite can be found \ref faq_pkcomposite "here
**/
namespace crule{
@@ -297,7 +298,7 @@ int main(int argc, char *argv[])
std::cout << "Error: resampling method " << resample_opt[0] << " not supported" << std::endl;
exit(1);
}
-
+
if(input_opt.empty()){
std::cerr << "No input file provided (use option -i). Use --help for help information" << std::endl;
exit(0);
@@ -306,7 +307,7 @@ int main(int argc, char *argv[])
int nwriteBand=0;
int writeBand=0;
vector<short> bands;
-
+
//get bounding box
double maxLRX=lrx_opt[0];
double maxULY=uly_opt[0];
@@ -349,20 +350,20 @@ int main(int argc, char *argv[])
exit(1);
}
if(!iextent){
- ulx_opt[0]=e_ulx;
- uly_opt[0]=e_uly;
- lrx_opt[0]=e_lrx;
- lry_opt[0]=e_lry;
+ ulx_opt[0]=e_ulx;
+ uly_opt[0]=e_uly;
+ lrx_opt[0]=e_lrx;
+ lry_opt[0]=e_lry;
}
else{
- if(e_ulx<ulx_opt[0])
- ulx_opt[0]=e_ulx;
- if(e_uly>uly_opt[0])
- uly_opt[0]=e_uly;
- if(e_lrx>lrx_opt[0])
- lrx_opt[0]=e_lrx;
- if(e_lry<lry_opt[0])
- lry_opt[0]=e_lry;
+ if(e_ulx<ulx_opt[0])
+ ulx_opt[0]=e_ulx;
+ if(e_uly>uly_opt[0])
+ uly_opt[0]=e_uly;
+ if(e_lrx>lrx_opt[0])
+ lrx_opt[0]=e_lrx;
+ if(e_lry<lry_opt[0])
+ lry_opt[0]=e_lry;
}
extentReader.close();
}
@@ -402,7 +403,7 @@ int main(int argc, char *argv[])
//todo: must be in init part only?
if(colorTable_opt.empty())
if(imgReader[ifile].getColorTable())
- theColorTable=(imgReader[ifile].getColorTable()->Clone());
+ theColorTable=(imgReader[ifile].getColorTable()->Clone());
if(projection_opt.empty())
theProjection=imgReader[ifile].getProjection();
if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
@@ -413,7 +414,7 @@ int main(int argc, char *argv[])
if((ulx_opt[0]||uly_opt[0]||lrx_opt[0]||lry_opt[0])&&(!imgReader[ifile].covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
if(verbose_opt[0])
- cout << input_opt[ifile] << " not within bounding box, skipping..." << endl;
+ cout << input_opt[ifile] << " not within bounding box, skipping..." << endl;
continue;
}
double theULX, theULY, theLRX, theLRY;
@@ -430,7 +431,7 @@ int main(int argc, char *argv[])
default:
case(crule::overwrite):
cout << "Composite rule: overwrite" << endl;
- break;
+ break;
case(crule::maxndvi):
cout << "Composite rule: max ndvi" << endl;
break;
@@ -467,7 +468,7 @@ int main(int argc, char *argv[])
}
}
if(band_opt.size()){
- nband=band_opt.size();
+ nband=band_opt.size();
bands.resize(band_opt.size());
for(int iband=0;iband<band_opt.size();++iband){
bands[iband]=band_opt[iband];
@@ -475,7 +476,7 @@ int main(int argc, char *argv[])
}
}
else{
- nband=imgReader[ifile].nrOfBand();
+ nband=imgReader[ifile].nrOfBand();
bands.resize(nband);
for(int iband=0;iband<nband;++iband)
bands[iband]=iband;
@@ -499,17 +500,17 @@ int main(int argc, char *argv[])
cout << "type of data for " << input_opt[ifile] << ": " << theType << endl;
cout << "nband: " << nband << endl;
}
-
+
maxLRX=theLRX;
maxULY=theULY;
minULX=theULX;
minLRY=theLRY;
if(dx_opt.size())
- dx=dx_opt[0];
+ dx=dx_opt[0];
else
dx=imgReader[ifile].getDeltaX();
if(dy_opt.size())
- dy=dy_opt[0];
+ dy=dy_opt[0];
else
dy=imgReader[ifile].getDeltaY();
init=true;
@@ -530,7 +531,7 @@ int main(int argc, char *argv[])
minULX=ulx_opt[0];
minLRY=lry_opt[0];
}
-
+
bool forceEUgrid=false;
if(projection_opt.size())
forceEUgrid=(!(projection_opt[0].compare("EPSG:3035"))||!(projection_opt[0].compare("EPSG:3035"))||projection_opt[0].find("ETRS-LAEA")!=string::npos);
@@ -573,14 +574,15 @@ int main(int argc, char *argv[])
//initialize image
if(verbose_opt[0])
cout << "initializing composite image..." << endl;
-// double dcol=(maxLRX-minULX+dx-1)/dx;
-// double drow=(maxULY-minLRY+dy-1)/dy;
-// int ncol=static_cast<int>(dcol);
-// int nrow=static_cast<int>(drow);
-
- int ncol=ceil((maxLRX-minULX)/dx);
- int nrow=ceil((maxULY-minLRY)/dy);
+ // double dcol=(maxLRX-minULX+dx-1)/dx;
+ // double drow=(maxULY-minLRY+dy-1)/dy;
+ // int ncol=static_cast<int>(dcol);
+ // int nrow=static_cast<int>(drow);
+ int ncol=abs(static_cast<int>(ceil((maxLRX-minULX)/dx)));
+ int nrow=abs(static_cast<int>(ceil((maxULY-minLRY)/dy)));
+ // int ncol=ceil((maxLRX-minULX)/dx);
+ // int nrow=ceil((maxULY-minLRY)/dy);
if(verbose_opt[0])
cout << "composite image dim (nrow x ncol): " << nrow << " x " << ncol << endl;
ImgWriterGdal imgWriter;
@@ -632,7 +634,7 @@ int main(int argc, char *argv[])
if(imgWriter.getDataType()==GDT_Byte){
if(colorTable_opt.size()){
if(colorTable_opt[0]!="none")
- imgWriter.setColorTable(colorTable_opt[0]);
+ imgWriter.setColorTable(colorTable_opt[0]);
}
else if(theColorTable)
imgWriter.setColorTable(theColorTable);
@@ -651,11 +653,11 @@ int main(int argc, char *argv[])
gt[5]=-dy;
maskWriter.setGeoTransform(gt);
if(projection_opt.size())
- maskWriter.setProjectionProj4(projection_opt[0]);
+ maskWriter.setProjectionProj4(projection_opt[0]);
else if(theProjection!=""){
- if(verbose_opt[0])
- cout << "projection: " << theProjection << endl;
- maskWriter.setProjection(theProjection);
+ if(verbose_opt[0])
+ cout << "projection: " << theProjection << endl;
+ maskWriter.setProjection(theProjection);
}
vector<double> burnValues(1,1);//burn value is 1 (single band)
maskWriter.rasterizeOgr(extentReader,burnValues,eoption_opt);
@@ -677,11 +679,11 @@ int main(int argc, char *argv[])
if(mask_opt.size()){
try{
if(verbose_opt[0]>=1)
- std::cout << "opening mask image file " << mask_opt[0] << std::endl;
+ std::cout << "opening mask image file " << mask_opt[0] << std::endl;
maskReader.open(mask_opt[0],GA_ReadOnly);
if(mskband_opt[0]>=maskReader.nrOfBand()){
- string errorString="Error: illegal mask band";
- throw(errorString);
+ string errorString="Error: illegal mask band";
+ throw(errorString);
}
}
catch(string error){
@@ -696,12 +698,12 @@ int main(int argc, char *argv[])
//create composite image
if(verbose_opt[0])
- cout << "creating composite image" << endl;
+ cout << "creating composite image" << endl;
Vector2d<double> writeBuffer(nband,imgWriter.nrOfCol());
vector<short> fileBuffer(ncol);//holds the number of used files
Vector2d<short> maxBuffer;//buffer used for maximum voting
// Vector2d<double> readBuffer(nband);
- vector<Vector2d<unsigned short> > readBuffer(input_opt.size());
+ vector<Vector2d<double> > readBuffer(input_opt.size());
for(int ifile=0;ifile<input_opt.size();++ifile)
readBuffer[ifile].resize(imgReader[ifile].nrOfBand());
statfactory::StatFactory stat;
@@ -776,7 +778,6 @@ int main(int argc, char *argv[])
ulj=floor(ulj);
lri=floor(lri);
lrj=floor(lrj);
-
double startCol=uli;
double endCol=lri;
if(uli<0)
@@ -787,7 +788,7 @@ int main(int argc, char *argv[])
endCol=0;
else if(lri>=imgReader[ifile].nrOfCol())
endCol=imgReader[ifile].nrOfCol()-1;
- int readncol=endCol-startCol+1;
+ // int readncol=endCol-startCol+1;
//lookup corresponding row for irow in this file
imgReader[ifile].geo2image(x,y,readCol,readRow);
@@ -797,59 +798,59 @@ int main(int argc, char *argv[])
}
// for(int iband=0;iband<imgReader.nrOfBand();++iband){
for(int iband=0;iband<nband;++iband){
- int readBand=(band_opt.size()>iband)? band_opt[iband] : iband;
+ int readBand=(band_opt.size()>iband)? band_opt[iband] : iband;
// readBuffer[iband].resize(readncol);
- try{
+ try{
imgReader[ifile].readData(readBuffer[ifile][iband],startCol,endCol,readRow,readBand,theResample);
- //test
+ //test
// imgReader[ifile].readData(readBuffer[ifile][iband],static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(readRow),static_cast<int>(readBand));
- // if(readRow==0&&iband==0){
- // for(int icol=0;icol<10;++icol)
- // cout << readBuffer[0][0][icol] << " ";
- // cout << endl;
- // }
- }
- catch(string error){
- cerr << "error reading image " << input_opt[ifile] << ": " << endl;
- throw;
- }
+ // if(readRow==0&&iband==0){
+ // for(int icol=0;icol<10;++icol)
+ // cout << readBuffer[0][0][icol] << " ";
+ // cout << endl;
+ // }
+ }
+ catch(string error){
+ cerr << "error reading image " << input_opt[ifile] << ": " << endl;
+ throw;
+ }
}
for(int ib=0;ib<ncol;++ib){
imgWriter.image2geo(ib,irow,x,y);
- //check mask first
- bool valid=true;
- if(mask_opt.size()){
- //read mask
- double colMask=0;
- double rowMask=0;
+ //check mask first
+ bool valid=true;
+ if(mask_opt.size()){
+ //read mask
+ double colMask=0;
+ double rowMask=0;
- maskReader.geo2image(x,y,colMask,rowMask);
- colMask=static_cast<int>(colMask);
- rowMask=static_cast<int>(rowMask);
- if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
- if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
+ maskReader.geo2image(x,y,colMask,rowMask);
+ colMask=static_cast<int>(colMask);
+ rowMask=static_cast<int>(rowMask);
+ if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
- assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
- try{
- maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
- }
- catch(string errorstring){
- cerr << errorstring << endl;
- exit(1);
- }
- catch(...){
- cerr << "error caught" << std::endl;
- exit(3);
- }
- oldRowMask=rowMask;
- }
- if(lineMask[colMask]==msknodata_opt[0])
- valid=false;
- }
- }
+ assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
+ try{
+ maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
+ }
+ catch(string errorstring){
+ cerr << errorstring << endl;
+ exit(1);
+ }
+ catch(...){
+ cerr << "error caught" << std::endl;
+ exit(3);
+ }
+ oldRowMask=rowMask;
+ }
+ if(lineMask[colMask]==msknodata_opt[0])
+ valid=false;
+ }
+ }
- if(!valid)
- continue;
+ if(!valid)
+ continue;
//lookup corresponding row for irow in this file
imgReader[ifile].geo2image(x,y,readCol,readRow);
@@ -870,64 +871,64 @@ int main(int argc, char *argv[])
upperCol=imgReader[ifile].nrOfCol()-1;
for(int vband=0;vband<bndnodata_opt.size();++vband){
val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][bndnodata_opt[vband]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][bndnodata_opt[vband]][lowerCol-startCol];
- if(minValue_opt.size()>vband){
- if(val_new<=minValue_opt[vband]){
- readValid=false;
- break;
- }
- }
- if(maxValue_opt.size()>vband){
- if(val_new>=maxValue_opt[vband]){
- readValid=false;
- break;
- }
- }
- if(srcnodata_opt.size()>vband){
- if(val_new==srcnodata_opt[vband]){
- readValid=false;
- break;
- }
- }
- }
+ if(minValue_opt.size()>vband){
+ if(val_new<=minValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(maxValue_opt.size()>vband){
+ if(val_new>=maxValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(srcnodata_opt.size()>vband){
+ if(val_new==srcnodata_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ }
break;
default:
readCol=static_cast<int>(readCol);
for(int vband=0;vband<bndnodata_opt.size();++vband){
val_new=readBuffer[ifile][bndnodata_opt[vband]][readCol-startCol];
- if(minValue_opt.size()>vband){
- if(val_new<=minValue_opt[vband]){
- readValid=false;
- break;
- }
- }
- if(maxValue_opt.size()>vband){
- if(val_new>=maxValue_opt[vband]){
- readValid=false;
- break;
- }
- }
- if(srcnodata_opt.size()>vband){
- if(val_new==srcnodata_opt[vband]){
- readValid=false;
- break;
- }
- }
- }
+ if(minValue_opt.size()>vband){
+ if(val_new<=minValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(maxValue_opt.size()>vband){
+ if(val_new>=maxValue_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ if(srcnodata_opt.size()>vband){
+ if(val_new==srcnodata_opt[vband]){
+ readValid=false;
+ break;
+ }
+ }
+ }
break;
- }
- if(readValid){
- if(file_opt[0]==1)
- ++fileBuffer[ib];
+ }
+ if(readValid){
+ if(file_opt[0]==1)
+ ++fileBuffer[ib];
if(writeValid[ib]){
int iband=0;
- switch(cruleMap[crule_opt[0]]){
- case(crule::maxndvi):{//max ndvi
+ switch(cruleMap[crule_opt[0]]){
+ case(crule::maxndvi):{//max ndvi
double red_current=writeBuffer[ruleBand_opt[0]][ib];
double nir_current=writeBuffer[ruleBand_opt[1]][ib];
- double ndvi_current=0;
+ double ndvi_current=0;
if(red_current+nir_current>0&&red_current>=0&&nir_current>=0)
ndvi_current=(nir_current-red_current)/(nir_current+red_current);
- double ndvi_new=0;
+ double ndvi_new=0;
double red_new=0;
double nir_new=0;
switch(theResample){
@@ -949,8 +950,8 @@ int main(int argc, char *argv[])
val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
writeBuffer[iband][ib]=val_new;
}
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
}
break;
default:
@@ -964,55 +965,55 @@ int main(int argc, char *argv[])
val_new=readBuffer[ifile][iband][readCol-startCol];
writeBuffer[iband][ib]=val_new;
}
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
}
break;
}
- break;
+ break;
}
- case(crule::maxband):
+ case(crule::maxband):
case(crule::minband):
case(crule::validband)://max,min,valid band
val_current=writeBuffer[ruleBand_opt[0]][ib];
- switch(theResample){
- case(BILINEAR):
- lowerCol=readCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=readCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReader[ifile].nrOfCol())
- upperCol=imgReader[ifile].nrOfCol()-1;
- val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][ruleBand_opt[0]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][ruleBand_opt[0]][lowerCol-startCol];
- val_new*=weight_opt[ifile];
- if((cruleMap[crule_opt[0]]==crule::maxband&&val_new>val_current)||(cruleMap[crule_opt[0]]==crule::minband&&val_new<val_current)||(cruleMap[crule_opt[0]]==crule::validband)){//&&val_new>minValue_opt[0]&&val_new<maxValue_opt[0])){
- for(iband=0;iband<nband;++iband){
- val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
- val_new*=weight_opt[ifile];
- writeBuffer[iband][ib]=val_new;
- }
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
+ switch(theResample){
+ case(BILINEAR):
+ lowerCol=readCol-0.5;
+ lowerCol=static_cast<int>(lowerCol);
+ upperCol=readCol+0.5;
+ upperCol=static_cast<int>(upperCol);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(upperCol>=imgReader[ifile].nrOfCol())
+ upperCol=imgReader[ifile].nrOfCol()-1;
+ val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][ruleBand_opt[0]][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][ruleBand_opt[0]][lowerCol-startCol];
+ val_new*=weight_opt[ifile];
+ if((cruleMap[crule_opt[0]]==crule::maxband&&val_new>val_current)||(cruleMap[crule_opt[0]]==crule::minband&&val_new<val_current)||(cruleMap[crule_opt[0]]==crule::validband)){//&&val_new>minValue_opt[0]&&val_new<maxValue_opt[0])){
+ for(iband=0;iband<nband;++iband){
+ val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
+ val_new*=weight_opt[ifile];
+ writeBuffer[iband][ib]=val_new;
}
- break;
- default:
- readCol=static_cast<int>(readCol);
- val_new=readBuffer[ifile][ruleBand_opt[0]][readCol-startCol];
- val_new*=weight_opt[ifile];
- if((cruleMap[crule_opt[0]]==crule::maxband&&val_new>val_current)||(cruleMap[crule_opt[0]]==crule::minband&&val_new<val_current)||(cruleMap[crule_opt[0]]==crule::validband)){//&&val_new>minValue_opt[0]&&val_new<maxValue_opt[0])){
- for(iband=0;iband<nband;++iband){
- val_new=readBuffer[ifile][iband][readCol-startCol];
- val_new*=weight_opt[ifile];
- writeBuffer[iband][ib]=val_new;
- }
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
+ }
+ break;
+ default:
+ readCol=static_cast<int>(readCol);
+ val_new=readBuffer[ifile][ruleBand_opt[0]][readCol-startCol];
+ val_new*=weight_opt[ifile];
+ if((cruleMap[crule_opt[0]]==crule::maxband&&val_new>val_current)||(cruleMap[crule_opt[0]]==crule::minband&&val_new<val_current)||(cruleMap[crule_opt[0]]==crule::validband)){//&&val_new>minValue_opt[0]&&val_new<maxValue_opt[0])){
+ for(iband=0;iband<nband;++iband){
+ val_new=readBuffer[ifile][iband][readCol-startCol];
+ val_new*=weight_opt[ifile];
+ writeBuffer[iband][ib]=val_new;
}
- break;
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
}
- break;
+ break;
+ }
+ break;
case(crule::mode)://max voting (only for Byte images)
switch(theResample){
case(BILINEAR):
@@ -1026,7 +1027,7 @@ int main(int argc, char *argv[])
upperCol=imgReader[ifile].nrOfCol()-1;
for(iband=0;iband<nband;++iband){
val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
- maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
+ maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
// ++(maxBuffer[ib][val_new]);
}
break;
@@ -1034,17 +1035,17 @@ int main(int argc, char *argv[])
readCol=static_cast<int>(readCol);
for(iband=0;iband<nband;++iband){
val_new=readBuffer[ifile][iband][readCol-startCol];
- maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
- }
+ maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
+ }
break;
- }
+ }
break;
case(crule::mean)://mean value
- case(crule::median)://median value
- case(crule::sum)://sum value
- case(crule::minallbands)://minimum for each and every band
- case(crule::maxallbands)://maximum for each and every band
- case(crule::stdev)://maximum for each and every band
+ case(crule::median)://median value
+ case(crule::sum)://sum value
+ case(crule::minallbands)://minimum for each and every band
+ case(crule::maxallbands)://maximum for each and every band
+ case(crule::stdev)://maximum for each and every band
switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
@@ -1073,11 +1074,11 @@ int main(int argc, char *argv[])
}
break;
}
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
- break;
- case(crule::overwrite):
- default:
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
+ break;
+ case(crule::overwrite):
+ default:
switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
@@ -1103,15 +1104,20 @@ int main(int argc, char *argv[])
}
break;
}
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
- break;
- }
- }
- else{
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
+ break;
+ }
+ }
+ else{
+ // //test
+ // if(readCol-startCol==0)
+ // std::cout << "val_new(" << ib << "): " << val_new << std::endl;
+ // else if(readCol-startCol==ncol-1)
+ // std::cout << "val_new(" << ib << "): " << val_new << std::endl;
writeValid[ib]=true;//readValid was true
int iband=0;
- switch(cruleMap[crule_opt[0]]){
+ switch(cruleMap[crule_opt[0]]){
case(crule::mean):
case(crule::median):
case(crule::sum):
@@ -1143,9 +1149,9 @@ int main(int argc, char *argv[])
}
break;
}
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
- break;
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
+ break;
case(crule::mode):
switch(theResample){
case(BILINEAR):
@@ -1159,17 +1165,17 @@ int main(int argc, char *argv[])
upperCol=imgReader[ifile].nrOfCol()-1;
for(iband=0;iband<nband;++iband){
val_new=(readCol-0.5-lowerCol)*readBuffer[ifile][iband][upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[ifile][iband][lowerCol-startCol];
- maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
+ maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
// ++(maxBuffer[ib][val_new]);
- }
+ }
break;
default:
readCol=static_cast<int>(readCol);
for(iband=0;iband<nband;++iband){
- val_new=readBuffer[ifile][iband][readCol-startCol];
- maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
- }
- // ++(maxBuffer[ib][val_new]);
+ val_new=readBuffer[ifile][iband][readCol-startCol];
+ maxBuffer[ib][val_new]=maxBuffer[ib][val_new]+weight_opt[ifile];
+ }
+ // ++(maxBuffer[ib][val_new]);
break;
}
break;
@@ -1199,8 +1205,8 @@ int main(int argc, char *argv[])
}
break;
}
- if(file_opt[0]>1)
- fileBuffer[ib]=ifile;
+ if(file_opt[0]>1)
+ fileBuffer[ib]=ifile;
break;
}
}
@@ -1228,8 +1234,8 @@ int main(int argc, char *argv[])
vector<short>::iterator maxit=maxBuffer[icol].begin();
maxit=stat.mymax(maxBuffer[icol],maxBuffer[icol].begin(),maxBuffer[icol].end());
writeBuffer[0][icol]=distance(maxBuffer[icol].begin(),maxit);
- if(file_opt[0]>1)
- fileBuffer[icol]=*(maxit);
+ if(file_opt[0]>1)
+ fileBuffer[icol]=*(maxit);
}
try{
imgWriter.writeData(writeBuffer[0],irow,0);
@@ -1247,42 +1253,42 @@ int main(int argc, char *argv[])
// assert(writeBuffer[bands[iband]].size()==imgWriter.nrOfCol());
assert(writeBuffer[iband].size()==imgWriter.nrOfCol());
for(int icol=0;icol<imgWriter.nrOfCol();++icol){
- try{
- switch(cruleMap[crule_opt[0]]){
- case(crule::mean):
- // writeBuffer[iband][icol]=stat.mean(storeBuffer[bands[iband]][icol]);
- writeBuffer[iband][icol]=stat.mean(storeBuffer[iband][icol]);
- break;
- case(crule::median):
- // writeBuffer[iband][icol]=stat.median(storeBuffer[bands[iband]][icol]);
- writeBuffer[iband][icol]=stat.median(storeBuffer[iband][icol]);
- break;
- case(crule::sum):
- // writeBuffer[iband][icol]=stat.sum(storeBuffer[bands[iband]][icol]);
- writeBuffer[iband][icol]=stat.sum(storeBuffer[iband][icol]);
- break;
- case(crule::minallbands):
- // writeBuffer[iband][icol]=stat.mymin(storeBuffer[bands[iband]][icol]);
- writeBuffer[iband][icol]=stat.mymin(storeBuffer[iband][icol]);
- break;
- case(crule::maxallbands):
- // writeBuffer[iband][icol]=stat.mymax(storeBuffer[bands[iband]][icol]);
- writeBuffer[iband][icol]=stat.mymax(storeBuffer[iband][icol]);
- break;
- case(crule::stdev):
- // writeBuffer[iband][icol]=sqrt(stat.var(storeBuffer[bands[iband]][icol]));
- writeBuffer[iband][icol]=sqrt(stat.var(storeBuffer[iband][icol]));
- break;
- default:
- break;
- }
- }
- catch(string error){
- if(verbose_opt[0])
- cerr << error << endl;
- writeBuffer[iband][icol]=dstnodata_opt[0];
- continue;
- }
+ try{
+ switch(cruleMap[crule_opt[0]]){
+ case(crule::mean):
+ // writeBuffer[iband][icol]=stat.mean(storeBuffer[bands[iband]][icol]);
+ writeBuffer[iband][icol]=stat.mean(storeBuffer[iband][icol]);
+ break;
+ case(crule::median):
+ // writeBuffer[iband][icol]=stat.median(storeBuffer[bands[iband]][icol]);
+ writeBuffer[iband][icol]=stat.median(storeBuffer[iband][icol]);
+ break;
+ case(crule::sum):
+ // writeBuffer[iband][icol]=stat.sum(storeBuffer[bands[iband]][icol]);
+ writeBuffer[iband][icol]=stat.sum(storeBuffer[iband][icol]);
+ break;
+ case(crule::minallbands):
+ // writeBuffer[iband][icol]=stat.mymin(storeBuffer[bands[iband]][icol]);
+ writeBuffer[iband][icol]=stat.mymin(storeBuffer[iband][icol]);
+ break;
+ case(crule::maxallbands):
+ // writeBuffer[iband][icol]=stat.mymax(storeBuffer[bands[iband]][icol]);
+ writeBuffer[iband][icol]=stat.mymax(storeBuffer[iband][icol]);
+ break;
+ case(crule::stdev):
+ // writeBuffer[iband][icol]=sqrt(stat.var(storeBuffer[bands[iband]][icol]));
+ writeBuffer[iband][icol]=sqrt(stat.var(storeBuffer[iband][icol]));
+ break;
+ default:
+ break;
+ }
+ }
+ catch(string error){
+ if(verbose_opt[0])
+ cerr << error << endl;
+ writeBuffer[iband][icol]=dstnodata_opt[0];
+ continue;
+ }
}
try{
imgWriter.writeData(writeBuffer[iband],irow,iband);
@@ -1315,4 +1321,3 @@ int main(int argc, char *argv[])
maskReader.close();
imgWriter.close();
}
-
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index f1d3d35..a0794bb 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -16,7 +16,7 @@ GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pktools. If not, see <http://www.gnu.org/licenses/>.
- ***********************************************************************/
+***********************************************************************/
#include <assert.h>
#include <cstdlib>
#include <string>
@@ -31,69 +31,69 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
/******************************************************************************/
/*! \page pkcrop pkcrop
- perform raster data operations on image such as crop, extract and stack bands
-## SYNOPSIS
+ perform raster data operations on image such as crop, extract and stack bands
+ ## SYNOPSIS
-<code>
+ <code>
Usage: pkcrop -i input -o output
-</code>
+ </code>
-<code>
+ <code>
- Options: [-of out_format] [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}] [-b band]* [-ulx ULX -uly ULY -lrx LRX -lry LRY] [-dx xres] [-dy yres] [-r resampling_method] [-a_srs epsg:number] [-nodata value]
+ Options: [-of out_format] [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}] [-b band]* [-ulx ULX -uly ULY -lrx LRX -lry LRY] [-dx xres] [-dy yres] [-r resampling_method] [-a_srs epsg:number] [-nodata value]
Advanced options:
- [-e vector [-cut] [-eo option]*] [-sband band -eband band]* [-co NAME=VALUE]* [-x center_x -y center_y] [-nx size_x -ny size_y] [-ns nsample -nl nlines] [-as min -as max] [-scale value]* [-off offset]* [-ct colortable] [-d description] [-align]
-</code>
-
-\section pkcrop_description Description
-
-The utility pkcrop can subset and stack raster images. In the spatial domain it can crop a bounding box from a larger image. The output bounding box is selected by setting the new corner coordinates using the options -ulx -uly -lrx -lry. Alternatively you can set the new image center (-x -y) and size. This can be done either in projected coordinates (using the options -nx -ny) or in image coordinates (using the options -ns -nl). You can also use a vector file to set the new bounding box [...]
-
-\section pkcrop_options Options
- - use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
- - short option `-h` shows basic options only, long option `--help` shows all options
-|short|long|type|default|description|
-|-----|----|----|-------|-----------|
- | i | input | std::string | |Input image file(s). If input contains multiple images, a multi-band output is created |
- | o | output | std::string | |Output image file |
- | a_srs | a_srs | std::string | |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 |
- | ulx | ulx | double | 0 |Upper left x value bounding box |
- | uly | uly | double | 0 |Upper left y value bounding box |
- | lrx | lrx | double | 0 |Lower right x value bounding box |
- | lry | lry | double | 0 |Lower right y value bounding box |
- | b | band | unsigned short | |band index to crop (leave empty to retain all bands) |
- | sband | startband | unsigned short | |Start band sequence number |
- | eband | endband | unsigned short | |End band sequence number |
- | as | autoscale | double | |scale output to min and max, e.g., --autoscale 0 --autoscale 255 |
- | ot | otype | std::string | |Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image |
- | of | oformat | std::string | GTiff |Output image format (see also gdal_translate)|
- | ct | ct | std::string | |color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid) |
- | dx | dx | double | |Output resolution in x (in meter) (empty: keep original resolution) |
- | dy | dy | double | |Output resolution in y (in meter) (empty: keep original resolution) |
- | r | resampling-method | std::string | near |Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation). |
- | e | extent | std::string | |get boundary from extent from polygons in vector file |
- | cut | crop_to_cutline | bool | false |Crop the extent of the target dataset to the extent of the cutline |
- | eo | eo | std::string | |special extent options controlling rasterization: ATTRIBUTE or CHUNKYSIZE or ALL_TOUCHED or BURN_VALUE_FROM or MERGE_ALG, e.g., -eo ATTRIBUTE=fieldname |
- | m | mask | std::string | |Use the specified file as a validity mask (0 is nodata) |
- | msknodata | msknodata | float | 0 |Mask value not to consider for crop
- | mskband | mskband | short | 0 |Mask band to read (0 indexed) |
- | co | co | std::string | |Creation option for output file. Multiple options can be specified. |
- | x | x | double | |x-coordinate of image center to crop (in meter) |
- | y | y | double | |y-coordinate of image center to crop (in meter) |
- | nx | nx | double | |image size in x to crop (in meter) |
- | ny | ny | double | |image size in y to crop (in meter) |
- | ns | ns | int | |number of samples to crop (in pixels) |
- | nl | nl | int | |number of lines to crop (in pixels) |
- | scale | scale | double | |output=scale*input+offset |
- | off | offset | double | |output=scale*input+offset |
- | nodata | nodata | float | |Nodata value to put in image if out of bounds. |
- | align | align | bool | |Align output bounding box to input image |
- | d | description | std::string | |Set image description |
-
-Examples
-========
-Some examples how to use pkcrop can be found \ref examples_pkcrop "here"
+ [-e vector [-cut] [-eo option]*] [-sband band -eband band]* [-co NAME=VALUE]* [-x center_x -y center_y] [-nx size_x -ny size_y] [-ns nsample -nl nlines] [-as min -as max] [-scale value]* [-off offset]* [-ct colortable] [-d description] [-align]
+ </code>
+
+ \section pkcrop_description Description
+
+ The utility pkcrop can subset and stack raster images. In the spatial domain it can crop a bounding box from a larger image. The output bounding box is selected by setting the new corner coordinates using the options -ulx -uly -lrx -lry. Alternatively you can set the new image center (-x -y) and size. This can be done either in projected coordinates (using the options -nx -ny) or in image coordinates (using the options -ns -nl). You can also use a vector file to set the new bounding bo [...]
+
+ \section pkcrop_options Options
+ - use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
+ - short option `-h` shows basic options only, long option `--help` shows all options
+ |short|long|type|default|description|
+ |-----|----|----|-------|-----------|
+ | i | input | std::string | |Input image file(s). If input contains multiple images, a multi-band output is created |
+ | o | output | std::string | |Output image file |
+ | a_srs | a_srs | std::string | |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 |
+ | ulx | ulx | double | 0 |Upper left x value bounding box |
+ | uly | uly | double | 0 |Upper left y value bounding box |
+ | lrx | lrx | double | 0 |Lower right x value bounding box |
+ | lry | lry | double | 0 |Lower right y value bounding box |
+ | b | band | unsigned short | |band index to crop (leave empty to retain all bands) |
+ | sband | startband | unsigned short | |Start band sequence number |
+ | eband | endband | unsigned short | |End band sequence number |
+ | as | autoscale | double | |scale output to min and max, e.g., --autoscale 0 --autoscale 255 |
+ | ot | otype | std::string | |Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image |
+ | of | oformat | std::string | GTiff |Output image format (see also gdal_translate)|
+ | ct | ct | std::string | |color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid) |
+ | dx | dx | double | |Output resolution in x (in meter) (empty: keep original resolution) |
+ | dy | dy | double | |Output resolution in y (in meter) (empty: keep original resolution) |
+ | r | resampling-method | std::string | near |Resampling method (near: nearest neighbor, bilinear: bi-linear interpolation). |
+ | e | extent | std::string | |get boundary from extent from polygons in vector file |
+ | cut | crop_to_cutline | bool | false |Crop the extent of the target dataset to the extent of the cutline |
+ | eo | eo | std::string | |special extent options controlling rasterization: ATTRIBUTE or CHUNKYSIZE or ALL_TOUCHED or BURN_VALUE_FROM or MERGE_ALG, e.g., -eo ATTRIBUTE=fieldname |
+ | m | mask | std::string | |Use the specified file as a validity mask (0 is nodata) |
+ | msknodata | msknodata | float | 0 |Mask value not to consider for crop
+ | mskband | mskband | short | 0 |Mask band to read (0 indexed) |
+ | co | co | std::string | |Creation option for output file. Multiple options can be specified. |
+ | x | x | double | |x-coordinate of image center to crop (in meter) |
+ | y | y | double | |y-coordinate of image center to crop (in meter) |
+ | nx | nx | double | |image size in x to crop (in meter) |
+ | ny | ny | double | |image size in y to crop (in meter) |
+ | ns | ns | int | |number of samples to crop (in pixels) |
+ | nl | nl | int | |number of lines to crop (in pixels) |
+ | scale | scale | double | |output=scale*input+offset |
+ | off | offset | double | |output=scale*input+offset |
+ | nodata | nodata | float | |Nodata value to put in image if out of bounds. |
+ | align | align | bool | |Align output bounding box to input image |
+ | d | description | std::string | |Set image description |
+
+ Examples
+ ========
+ Some examples how to use pkcrop can be found \ref examples_pkcrop "here"
**/
using namespace std;
@@ -123,8 +123,8 @@ int main(int argc, char *argv[])
Optionpk<int> ns_opt("ns", "ns", "number of samples to crop (in pixels)");
Optionpk<int> nl_opt("nl", "nl", "number of lines to crop (in pixels)");
Optionpk<unsigned short> band_opt("b", "band", "band index to crop (leave empty to retain all bands)");
- Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
- Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
+ Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band sequence number");
+ Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence number");
Optionpk<double> autoscale_opt("as", "autoscale", "scale output to min and max, e.g., --autoscale 0 --autoscale 255");
Optionpk<double> scale_opt("scale", "scale", "output=scale*input+offset");
Optionpk<double> offset_opt("offset", "offset", "output=scale*input+offset");
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
offset_opt.setHide(1);
nodata_opt.setHide(1);
description_opt.setHide(1);
-
+
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
@@ -257,17 +257,17 @@ int main(int argc, char *argv[])
try{
if(bstart_opt.size()){
if(bend_opt.size()!=bstart_opt.size()){
- string errorstring="Error: options for start and end band indexes must be provided as pairs, missing end band";
- throw(errorstring);
+ string errorstring="Error: options for start and end band indexes must be provided as pairs, missing end band";
+ throw(errorstring);
}
band_opt.clear();
for(int ipair=0;ipair<bstart_opt.size();++ipair){
- if(bend_opt[ipair]<=bstart_opt[ipair]){
- string errorstring="Error: index for end band must be smaller then start band";
- throw(errorstring);
- }
- for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
- band_opt.push_back(iband);
+ if(bend_opt[ipair]<=bstart_opt[ipair]){
+ string errorstring="Error: index for end band must be smaller then start band";
+ throw(errorstring);
+ }
+ for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
+ band_opt.push_back(iband);
}
}
}
@@ -280,7 +280,7 @@ int main(int argc, char *argv[])
string projectionString;
for(int iimg=0;iimg<input_opt.size();++iimg){
try{
- imgReader.open(input_opt[iimg],GA_ReadOnly);
+ imgReader.open(input_opt[iimg],GA_ReadOnly);
}
catch(string error){
cerr << "Error: could not open file " << input_opt[iimg] << ": " << error << std::endl;
@@ -294,7 +294,7 @@ int main(int argc, char *argv[])
if(!iimg||imgReader.getDeltaX()<dx)
dx=imgReader.getDeltaX();
}
-
+
if(dy_opt.empty()){
if(!iimg||imgReader.getDeltaY()<dy)
dy=imgReader.getDeltaY();
@@ -351,20 +351,20 @@ int main(int argc, char *argv[])
exit(1);
}
if(!iextent){
- ulx_opt[0]=e_ulx;
- uly_opt[0]=e_uly;
- lrx_opt[0]=e_lrx;
- lry_opt[0]=e_lry;
+ ulx_opt[0]=e_ulx;
+ uly_opt[0]=e_uly;
+ lrx_opt[0]=e_lrx;
+ lry_opt[0]=e_lry;
}
else{
- if(e_ulx<ulx_opt[0])
- ulx_opt[0]=e_ulx;
- if(e_uly>uly_opt[0])
- uly_opt[0]=e_uly;
- if(e_lrx>lrx_opt[0])
- lrx_opt[0]=e_lrx;
- if(e_lry<lry_opt[0])
- lry_opt[0]=e_lry;
+ if(e_ulx<ulx_opt[0])
+ ulx_opt[0]=e_ulx;
+ if(e_uly>uly_opt[0])
+ uly_opt[0]=e_uly;
+ if(e_lrx>lrx_opt[0])
+ lrx_opt[0]=e_lrx;
+ if(e_lry<lry_opt[0])
+ lry_opt[0]=e_lry;
}
extentReader.close();
}
@@ -434,10 +434,10 @@ int main(int argc, char *argv[])
gt[5]=-dy;
maskWriter.setGeoTransform(gt);
if(projection_opt.size())
- maskWriter.setProjectionProj4(projection_opt[0]);
+ maskWriter.setProjectionProj4(projection_opt[0]);
else if(projectionString.size())
- maskWriter.setProjection(projectionString);
-
+ maskWriter.setProjection(projectionString);
+
vector<double> burnValues(1,1);//burn value is 1 (single band)
maskWriter.rasterizeOgr(extentReader,burnValues,eoption_opt);
maskWriter.close();
@@ -455,8 +455,8 @@ int main(int argc, char *argv[])
//there is only a single mask
maskReader.open(mask_opt[0],GA_ReadOnly);
if(mskband_opt[0]>=maskReader.nrOfBand()){
- string errorString="Error: illegal mask band";
- throw(errorString);
+ string errorString="Error: illegal mask band";
+ throw(errorString);
}
}
catch(string error){
@@ -528,12 +528,12 @@ int main(int argc, char *argv[])
double magicX=1,magicY=1;
// imgReader.getMagicPixel(magicX,magicY);
if(forceEUgrid){
- //force to LAEA grid
- Egcs egcs;
+ //force to LAEA grid
+ Egcs egcs;
egcs.setLevel(egcs.res2level(dx));
- egcs.force2grid(cropulx,cropuly,croplrx,croplry);
- imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
- imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
+ egcs.force2grid(cropulx,cropuly,croplrx,croplry);
+ imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
+ imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
}
imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
@@ -549,28 +549,28 @@ int main(int argc, char *argv[])
croplrx=lrx_opt[0];
croplry=lry_opt[0];
if(forceEUgrid){
- //force to LAEA grid
- Egcs egcs;
+ //force to LAEA grid
+ Egcs egcs;
egcs.setLevel(egcs.res2level(dx));
- egcs.force2grid(cropulx,cropuly,croplrx,croplry);
+ egcs.force2grid(cropulx,cropuly,croplrx,croplry);
}
else if(align_opt[0]){
- if(cropulx>imgReader.getUlx())
- cropulx-=fmod(cropulx-imgReader.getUlx(),dx);
- else if(cropulx<imgReader.getUlx())
- cropulx+=fmod(imgReader.getUlx()-cropulx,dx)-dx;
- if(croplrx<imgReader.getLrx())
- croplrx+=fmod(imgReader.getLrx()-croplrx,dx);
- else if(croplrx>imgReader.getLrx())
- croplrx-=fmod(croplrx-imgReader.getLrx(),dx)+dx;
- if(croplry>imgReader.getLry())
- croplry-=fmod(croplry-imgReader.getLry(),dy);
- else if(croplry<imgReader.getLry())
- croplry+=fmod(imgReader.getLry()-croplry,dy)-dy;
- if(cropuly<imgReader.getUly())
- cropuly+=fmod(imgReader.getUly()-cropuly,dy);
- else if(cropuly>imgReader.getUly())
- cropuly-=fmod(cropuly-imgReader.getUly(),dy)+dy;
+ if(cropulx>imgReader.getUlx())
+ cropulx-=fmod(cropulx-imgReader.getUlx(),dx);
+ else if(cropulx<imgReader.getUlx())
+ cropulx+=fmod(imgReader.getUlx()-cropulx,dx)-dx;
+ if(croplrx<imgReader.getLrx())
+ croplrx+=fmod(imgReader.getLrx()-croplrx,dx);
+ else if(croplrx>imgReader.getLrx())
+ croplrx-=fmod(croplrx-imgReader.getLrx(),dx)+dx;
+ if(croplry>imgReader.getLry())
+ croplry-=fmod(croplry-imgReader.getLry(),dy);
+ else if(croplry<imgReader.getLry())
+ croplry+=fmod(imgReader.getLry()-croplry,dy)-dy;
+ if(cropuly<imgReader.getUly())
+ cropuly+=fmod(imgReader.getUly()-cropuly,dy);
+ else if(cropuly>imgReader.getUly())
+ cropuly-=fmod(cropuly-imgReader.getUly(),dy)+dy;
}
imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
@@ -591,21 +591,21 @@ int main(int argc, char *argv[])
dcroprow=(lrj-ulj+1)/(dy/deltaY);
if(!imgWriter.nrOfBand()){//not opened yet
if(verbose_opt[0]){
- cout << "cropulx: " << cropulx << endl;
- cout << "cropuly: " << cropuly << endl;
- cout << "croplrx: " << croplrx << endl;
- cout << "croplry: " << croplry << endl;
- cout << "ncropcol: " << ncropcol << endl;
- cout << "ncroprow: " << ncroprow << endl;
- cout << "cropulx+ncropcol*dx: " << cropulx+ncropcol*dx << endl;
- cout << "cropuly-ncroprow*dy: " << cropuly-ncroprow*dy << endl;
- cout << "upper left column of input image: " << uli << endl;
- cout << "upper left row of input image: " << ulj << endl;
- cout << "lower right column of input image: " << lri << endl;
- cout << "lower right row of input image: " << lrj << endl;
- cout << "new number of cols: " << ncropcol << endl;
- cout << "new number of rows: " << ncroprow << endl;
- cout << "new number of bands: " << ncropband << endl;
+ cout << "cropulx: " << cropulx << endl;
+ cout << "cropuly: " << cropuly << endl;
+ cout << "croplrx: " << croplrx << endl;
+ cout << "croplry: " << croplry << endl;
+ cout << "ncropcol: " << ncropcol << endl;
+ cout << "ncroprow: " << ncroprow << endl;
+ cout << "cropulx+ncropcol*dx: " << cropulx+ncropcol*dx << endl;
+ cout << "cropuly-ncroprow*dy: " << cropuly-ncroprow*dy << endl;
+ cout << "upper left column of input image: " << uli << endl;
+ cout << "upper left row of input image: " << ulj << endl;
+ cout << "lower right column of input image: " << lri << endl;
+ cout << "lower right row of input image: " << lrj << endl;
+ cout << "new number of cols: " << ncropcol << endl;
+ cout << "new number of rows: " << ncroprow << endl;
+ cout << "new number of bands: " << ncropband << endl;
}
// string theCompression;
// if(compress_opt[0]!="")//default
@@ -622,18 +622,18 @@ int main(int argc, char *argv[])
imageType=oformat_opt[0];
try{
imgWriter.open(output_opt[0],ncropcol,ncroprow,ncropband,theType,imageType,option_opt);
- if(nodata_opt.size()){
- imgWriter.setNoData(nodata_opt);
- // for(int iband=0;iband<ncropband;++iband)
- // imgWriter.GDALSetNoDataValue(nodata_opt[0],iband);
- }
+ if(nodata_opt.size()){
+ imgWriter.setNoData(nodata_opt);
+ // for(int iband=0;iband<ncropband;++iband)
+ // imgWriter.GDALSetNoDataValue(nodata_opt[0],iband);
+ }
}
catch(string errorstring){
cout << errorstring << endl;
exit(4);
}
if(description_opt.size())
- imgWriter.setImageDescription(description_opt[0]);
+ imgWriter.setImageDescription(description_opt[0]);
double gt[6];
gt[0]=cropulx;
gt[1]=dx;
@@ -643,19 +643,19 @@ int main(int argc, char *argv[])
gt[5]=(imgReader.isGeoRef())? -dy : dy;
imgWriter.setGeoTransform(gt);
if(projection_opt.size()){
- if(verbose_opt[0])
- cout << "projection: " << projection_opt[0] << endl;
- imgWriter.setProjectionProj4(projection_opt[0]);
+ if(verbose_opt[0])
+ cout << "projection: " << projection_opt[0] << endl;
+ imgWriter.setProjectionProj4(projection_opt[0]);
}
else
- imgWriter.setProjection(imgReader.getProjection());
+ imgWriter.setProjection(imgReader.getProjection());
if(imgWriter.getDataType()==GDT_Byte){
- if(colorTable_opt.size()){
- if(colorTable_opt[0]!="none")
- imgWriter.setColorTable(colorTable_opt[0]);
- }
- else if (imgReader.getColorTable()!=NULL)//copy colorTable from input image
- imgWriter.setColorTable(imgReader.getColorTable());
+ if(colorTable_opt.size()){
+ if(colorTable_opt[0]!="none")
+ imgWriter.setColorTable(colorTable_opt[0]);
+ }
+ else if (imgReader.getColorTable()!=NULL)//copy colorTable from input image
+ imgWriter.setColorTable(imgReader.getColorTable());
}
}
@@ -680,8 +680,6 @@ int main(int argc, char *argv[])
else if(lrj>=imgReader.nrOfRow())
endRow=imgReader.nrOfRow()-1;
-
-
int readncol=endCol-startCol+1;
// vector<double> readBuffer(readncol+1);
vector<double> readBuffer;
@@ -689,38 +687,38 @@ int main(int argc, char *argv[])
for(int iband=0;iband<nband;++iband){
int readBand=(band_opt.size()>iband)?band_opt[iband]:iband;
if(verbose_opt[0]){
- cout << "extracting band " << readBand << endl;
- pfnProgress(progress,pszMessage,pProgressArg);
+ cout << "extracting band " << readBand << endl;
+ pfnProgress(progress,pszMessage,pProgressArg);
}
double theMin=0;
double theMax=0;
if(autoscale_opt.size()){
- try{
- imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax);
- }
- catch(string errorString){
- cout << errorString << endl;
- }
- if(verbose_opt[0])
- cout << "minmax: " << theMin << ", " << theMax << endl;
- double theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin);
- double theOffset=autoscale_opt[0]-theScale*theMin;
- imgReader.setScale(theScale,readBand);
- imgReader.setOffset(theOffset,readBand);
- }
+ try{
+ imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax);
+ }
+ catch(string errorString){
+ cout << errorString << endl;
+ }
+ if(verbose_opt[0])
+ cout << "minmax: " << theMin << ", " << theMax << endl;
+ double theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin);
+ double theOffset=autoscale_opt[0]-theScale*theMin;
+ imgReader.setScale(theScale,readBand);
+ imgReader.setOffset(theOffset,readBand);
+ }
else{
- if(scale_opt.size()){
- if(scale_opt.size()>iband)
- imgReader.setScale(scale_opt[iband],readBand);
- else
- imgReader.setScale(scale_opt[0],readBand);
- }
- if(offset_opt.size()){
- if(offset_opt.size()>iband)
- imgReader.setOffset(offset_opt[iband],readBand);
- else
- imgReader.setOffset(offset_opt[0],readBand);
- }
+ if(scale_opt.size()){
+ if(scale_opt.size()>iband)
+ imgReader.setScale(scale_opt[iband],readBand);
+ else
+ imgReader.setScale(scale_opt[0],readBand);
+ }
+ if(offset_opt.size()){
+ if(offset_opt.size()>iband)
+ imgReader.setOffset(offset_opt[iband],readBand);
+ else
+ imgReader.setOffset(offset_opt[0],readBand);
+ }
}
double readRow=0;
@@ -728,84 +726,84 @@ int main(int argc, char *argv[])
double lowerCol=0;
double upperCol=0;
for(int irow=0;irow<imgWriter.nrOfRow();++irow){
- vector<float> lineMask;
- double x=0;
- double y=0;
- //convert irow to geo
- imgWriter.image2geo(0,irow,x,y);
- //lookup corresponding row for irow in this file
- imgReader.geo2image(x,y,readCol,readRow);
- vector<double> writeBuffer;
- if(readRow<0||readRow>=imgReader.nrOfRow()){
- //if(readRow<0)
- //readRow=0;
- //else if(readRow>=imgReader.nrOfRow())
- //readRow=imgReader.nrOfRow()-1;
- for(int icol=0;icol<imgWriter.nrOfCol();++icol)
- writeBuffer.push_back(nodataValue);
- }
- else{
- try{
+ vector<float> lineMask;
+ double x=0;
+ double y=0;
+ //convert irow to geo
+ imgWriter.image2geo(0,irow,x,y);
+ //lookup corresponding row for irow in this file
+ imgReader.geo2image(x,y,readCol,readRow);
+ vector<double> writeBuffer;
+ if(readRow<0||readRow>=imgReader.nrOfRow()){
+ //if(readRow<0)
+ //readRow=0;
+ //else if(readRow>=imgReader.nrOfRow())
+ //readRow=imgReader.nrOfRow()-1;
+ for(int icol=0;icol<imgWriter.nrOfCol();++icol)
+ writeBuffer.push_back(nodataValue);
+ }
+ else{
+ try{
if(endCol<imgReader.nrOfCol()-1){
imgReader.readData(readBuffer,startCol,endCol+1,readRow,readBand,theResample);
}
else{
imgReader.readData(readBuffer,startCol,endCol,readRow,readBand,theResample);
}
- // for(int icol=0;icol<ncropcol;++icol){
- double oldRowMask=-1;//keep track of row mask to optimize number of line readings
- for(int icol=0;icol<imgWriter.nrOfCol();++icol){
- imgWriter.image2geo(icol,irow,x,y);
- //lookup corresponding row for irow in this file
- imgReader.geo2image(x,y,readCol,readRow);
- if(readCol<0||readCol>=imgReader.nrOfCol()){
- // if(readCol<0||readCol>=imgReader.nrOfCol()){
- // if(readCol<0)
- // readCol=0;
- // else if(readCol>=imgReader.nrOfCol())
- // readCol=imgReader.nrOfCol()-1;
- writeBuffer.push_back(nodataValue);
- }
- else{
-
+ // for(int icol=0;icol<ncropcol;++icol){
+ double oldRowMask=-1;//keep track of row mask to optimize number of line readings
+ for(int icol=0;icol<imgWriter.nrOfCol();++icol){
+ imgWriter.image2geo(icol,irow,x,y);
+ //lookup corresponding row for irow in this file
+ imgReader.geo2image(x,y,readCol,readRow);
+ if(readCol<0||readCol>=imgReader.nrOfCol()){
+ // if(readCol<0||readCol>=imgReader.nrOfCol()){
+ // if(readCol<0)
+ // readCol=0;
+ // else if(readCol>=imgReader.nrOfCol())
+ // readCol=imgReader.nrOfCol()-1;
+ writeBuffer.push_back(nodataValue);
+ }
+ else{
+
bool valid=true;
- double geox=0;
- double geoy=0;
+ double geox=0;
+ double geoy=0;
if(mask_opt.size()){
- //read mask
- double colMask=0;
- double rowMask=0;
-
- imgWriter.image2geo(icol,irow,geox,geoy);
- maskReader.geo2image(geox,geoy,colMask,rowMask);
- colMask=static_cast<int>(colMask);
- rowMask=static_cast<int>(rowMask);
- if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
- if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
-
- assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
- try{
- maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
- }
- catch(string errorstring){
- cerr << errorstring << endl;
- exit(1);
- }
- catch(...){
- cerr << "error caught" << std::endl;
- exit(3);
- }
- oldRowMask=rowMask;
- }
+ //read mask
+ double colMask=0;
+ double rowMask=0;
+
+ imgWriter.image2geo(icol,irow,geox,geoy);
+ maskReader.geo2image(geox,geoy,colMask,rowMask);
+ colMask=static_cast<int>(colMask);
+ rowMask=static_cast<int>(rowMask);
+ if(rowMask>=0&&rowMask<maskReader.nrOfRow()&&colMask>=0&&colMask<maskReader.nrOfCol()){
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldRowMask)){
+
+ assert(rowMask>=0&&rowMask<maskReader.nrOfRow());
+ try{
+ maskReader.readData(lineMask,static_cast<int>(rowMask),mskband_opt[0]);
+ }
+ catch(string errorstring){
+ cerr << errorstring << endl;
+ exit(1);
+ }
+ catch(...){
+ cerr << "error caught" << std::endl;
+ exit(3);
+ }
+ oldRowMask=rowMask;
+ }
for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
if(lineMask[colMask]==msknodata_opt[ivalue]){
valid=false;
- if(nodata_opt.size()>ivalue)
- nodataValue=nodata_opt[ivalue];
+ if(nodata_opt.size()>ivalue)
+ nodataValue=nodata_opt[ivalue];
}
}
- }
- }
+ }
+ }
if(!valid)
writeBuffer.push_back(nodataValue);
else{
@@ -828,39 +826,39 @@ int main(int argc, char *argv[])
// writeBuffer.push_back(readBuffer[readCol]*theScale+theOffset);
writeBuffer.push_back(readBuffer[readCol]);
break;
- }
- }
- }
- }
- }
- catch(string errorstring){
- cout << errorstring << endl;
- exit(2);
- }
- }
- if(writeBuffer.size()!=imgWriter.nrOfCol())
- cout << "writeBuffer.size()=" << writeBuffer.size() << ", imgWriter.nrOfCol()=" << imgWriter.nrOfCol() << endl;
- assert(writeBuffer.size()==imgWriter.nrOfCol());
- try{
- imgWriter.writeData(writeBuffer,irow,writeBand);
- }
- catch(string errorstring){
- cout << errorstring << endl;
- exit(3);
- }
- if(verbose_opt[0]){
- progress=(1.0+irow);
- progress/=imgWriter.nrOfRow();
- pfnProgress(progress,pszMessage,pProgressArg);
- }
- else{
- progress=(1.0+irow);
- progress+=(imgWriter.nrOfRow()*writeBand);
- progress/=imgWriter.nrOfBand()*imgWriter.nrOfRow();
- assert(progress>=0);
- assert(progress<=1);
- pfnProgress(progress,pszMessage,pProgressArg);
- }
+ }
+ }
+ }
+ }
+ }
+ catch(string errorstring){
+ cout << errorstring << endl;
+ exit(2);
+ }
+ }
+ if(writeBuffer.size()!=imgWriter.nrOfCol())
+ cout << "writeBuffer.size()=" << writeBuffer.size() << ", imgWriter.nrOfCol()=" << imgWriter.nrOfCol() << endl;
+ assert(writeBuffer.size()==imgWriter.nrOfCol());
+ try{
+ imgWriter.writeData(writeBuffer,irow,writeBand);
+ }
+ catch(string errorstring){
+ cout << errorstring << endl;
+ exit(3);
+ }
+ if(verbose_opt[0]){
+ progress=(1.0+irow);
+ progress/=imgWriter.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ else{
+ progress=(1.0+irow);
+ progress+=(imgWriter.nrOfRow()*writeBand);
+ progress/=imgWriter.nrOfBand()*imgWriter.nrOfRow();
+ assert(progress>=0);
+ assert(progress<=1);
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
}
++writeBand;
}
diff --git a/src/apps/pkkalman.cc b/src/apps/pkkalman.cc
index 344ffe2..2a111c8 100644
--- a/src/apps/pkkalman.cc
+++ b/src/apps/pkkalman.cc
@@ -45,43 +45,43 @@ produce kalman filtered raster time series
\section pkkalman_description Description
-The utilty pkkalman will complement a time series of observations (option -obs) at fine spatial resolution. A data assimilation technique based on a Kalman filter is hereby used. The data at fine spatial resolution are assimilated with coarse spatial resolution time series at a finer temporal resolution, referred to as a model (option -mod). The time series for both observation and model can either be provided as multi-band raster datasets or as multiple single band datasets. Missing dat [...]
+The utilty pkkalman will complement a time series of observations (option -obs) at fine spatial resolution. A data assimilation technique based on a Kalman filter is hereby used. The data at fine spatial resolution are assimilated with coarse spatial resolution time series at a finer temporal resolution, referred to as a model (option -mod). The time series for both observation and model can either be provided as multi-band raster datasets or as multiple single band datasets. Missing dat [...]
\section pkcrop_options Options
- use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
- short option `-h` shows basic options only, long option `--help` shows all options
|short|long|type|default|description|
|-----|----|----|-------|-----------|
- | dir | direction | std::string | forward |direction to run model (forward\|backward\|smooth) |
- | mod | model | std::string | |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.) |
- | modmask| modmask | std::string | |model mask datasets(s). Must have same dimension as model input. Use either multi-band input or multiple single-band inputs|
- | obs | observation | std::string | |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.) |
- | obsmask| obsmask | std::string | |observation mask dataset(s). Must have same dimension as observation input (use multi-band input or multiple single-band inputs |
- | tmod | tmodel | int | |time sequence of model input. Sequence must have exact same length as model input. Leave empty to have default sequence 0,1,2,etc. |
- | tobs | tobservation | int | |time sequence of observation input. Sequence must have exact same length as observation input |
- | a_srs | a_srs | std::string | |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 |
- | ofw | outputfw | std::string | |Output raster dataset for forward model |
- | u_ofw | u_outputfw | std::string | |Uncertainty output raster dataset for forward model |
- | obw | outputbw | std::string | |Output raster dataset for backward model |
- | u_obw | u_outputbw | std::string | |Uncertainty output raster dataset for backward model |
- | ofb | outputfb | std::string | |Output raster dataset for smooth model |
- | u_ofb | u_outputfb | std::string | |Uncertainty output raster dataset for smooth model |
- | modnodata | modnodata | double | 0 |invalid value for model input |
- | obsnodata | obsnodata | double | 0 |invalid value for observation input |
+ | dir | direction | std::string | forward |direction to run model (forward\|backward\|smooth) |
+ | mod | model | std::string | |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.) |
+ | modmask| modmask | std::string | |model mask datasets(s). Must have same dimension as model input. Use either multi-band input or multiple single-band inputs|
+ | obs | observation | std::string | |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.) |
+ | obsmask| obsmask | std::string | |observation mask dataset(s). Must have same dimension as observation input (use multi-band input or multiple single-band inputs |
+ | tmod | tmodel | int | |time sequence of model input. Sequence must have exact same length as model input. Leave empty to have default sequence 0,1,2,etc. |
+ | tobs | tobservation | int | |time sequence of observation input. Sequence must have exact same length as observation input |
+ | a_srs | a_srs | std::string | |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 |
+ | ofw | outputfw | std::string | |Output raster dataset for forward model |
+ | u_ofw | u_outputfw | std::string | |Uncertainty output raster dataset for forward model |
+ | obw | outputbw | std::string | |Output raster dataset for backward model |
+ | u_obw | u_outputbw | std::string | |Uncertainty output raster dataset for backward model |
+ | ofb | outputfb | std::string | |Output raster dataset for smooth model |
+ | u_ofb | u_outputfb | std::string | |Uncertainty output raster dataset for smooth model |
+ | modnodata | modnodata | double | 0 |invalid value for model input |
+ | obsnodata | obsnodata | double | 0 |invalid value for observation input |
| msknodata | msknodata | float | 0 |Mask value not to consider
- | mskband | mskband | short | 0 |Mask band to read (0 indexed) |
- | obsmin | obsmin | double | |Minimum value for observation data |
- | obsmax | obsmax | double | |Maximum value for observation data |
- | eps | eps | double | 1e-05 |epsilon for non zero division |
- | um | uncertmodel | double | 1 |Uncertainty of the model |
- | uo | uncertobs | double | 1 |Uncertainty of valid observations |
- | unodata | uncertnodata | double | 100 |Uncertainty in case of no-data values in observation |
- | q | q | double | 1 |Process noise: expresses instability (variance) of proportions of fine res pixels within a moderate resolution pixel |
- | down | down | int | |Downsampling factor for reading model data to calculate regression (default is ratio between coarse (model) and fine (obs) resolution raster datasets)|
- | ot | otype | std::string | |Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image |
- | of | oformat | std::string | GTiff |Output image format (see also gdal_translate).|
- | co | co | std::string | |Creation option for output file. Multiple options can be specified. |
- | v | verbose | short | 0 |verbose mode when positive |
+ | mskband | mskband | short | 0 |Mask band to read (0 indexed) |
+ | obsmin | obsmin | double | |Minimum value for observation data |
+ | obsmax | obsmax | double | |Maximum value for observation data |
+ | eps | eps | double | 1e-05 |epsilon for non zero division |
+ | um | uncertmodel | double | 1 |Uncertainty of the model |
+ | uo | uncertobs | double | 1 |Uncertainty of valid observations |
+ | unodata | uncertnodata | double | 100 |Uncertainty in case of no-data values in observation |
+ | q | q | double | 1 |Process noise: expresses instability (variance) of proportions of fine res pixels within a moderate resolution pixel |
+ | down | down | int | |Downsampling factor for reading model data to calculate regression (default is ratio between coarse (model) and fine (obs) resolution raster datasets)|
+ | ot | otype | std::string | |Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image |
+ | of | oformat | std::string | GTiff |Output image format (see also gdal_translate).|
+ | co | co | std::string | |Creation option for output file. Multiple options can be specified. |
+ | v | verbose | short | 0 |verbose mode when positive |
Examples
========
@@ -98,8 +98,8 @@ int main(int argc,char **argv) {
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");
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.)");
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");
- Optionpk<int> tmodel_opt("tmod","tmodel","time sequence of model input. Sequence must have exact same length as model input. Leave empty to have default sequence 0,1,2,etc.");
- Optionpk<int> tobservation_opt("tobs","tobservation","time sequence of observation input. Sequence must have exact same length as observation input)");
+ Optionpk<int> tmodel_opt("tmod","tmodel","time sequence of model input. Sequence must have exact same length as model input. Leave empty to have default sequence 0,1,2,etc.");
+ Optionpk<int> tobservation_opt("tobs","tobservation","time sequence of observation input. Sequence must have exact same length as observation input)");
Optionpk<string> projection_opt("a_srs", "a_srs", "Override the projection for the output file (leave blank to copy from input file, use epsg:3035 to use European projection and force to European grid");
Optionpk<string> outputfw_opt("ofw", "outputfw", "Output raster dataset for forward model");
Optionpk<string> uncertfw_opt("u_ofw", "u_outputfw", "Uncertainty output raster dataset for forward model");
@@ -201,29 +201,29 @@ int main(int argc,char **argv) {
}
if(find(direction_opt.begin(),direction_opt.end(),"forward")!=direction_opt.end()){
if(outputfw_opt.empty()){
- errorStream << "Error: output forward datasets is not provided, use option -ofw" << endl;
- throw(errorStream.str());
+ errorStream << "Error: output forward datasets is not provided, use option -ofw" << endl;
+ throw(errorStream.str());
}
if(uncertfw_opt.empty()){
- ostringstream uncertStream;
- uncertStream << outputfw_opt[0] << "_uncert";
- uncertfw_opt.push_back(uncertStream.str());
+ ostringstream uncertStream;
+ uncertStream << outputfw_opt[0] << "_uncert";
+ uncertfw_opt.push_back(uncertStream.str());
}
}
if(find(direction_opt.begin(),direction_opt.end(),"backward")!=direction_opt.end()){
if(outputbw_opt.empty()){
- errorStream << "Error: output backward datasets is not provided, use option -obw" << endl;
- throw(errorStream.str());
+ errorStream << "Error: output backward datasets is not provided, use option -obw" << endl;
+ throw(errorStream.str());
}
if(uncertbw_opt.empty()){
- ostringstream uncertStream;
- uncertStream << outputbw_opt[0] << "_uncert";
- uncertbw_opt.push_back(uncertStream.str());
+ ostringstream uncertStream;
+ uncertStream << outputbw_opt[0] << "_uncert";
+ uncertbw_opt.push_back(uncertStream.str());
}
}
// if(model_opt.size()<observation_opt.size()){
- // errorStream << "Error: sequence of models should be larger than observations" << endl;
- // throw(errorStream.str());
+ // errorStream << "Error: sequence of models should be larger than observations" << endl;
+ // throw(errorStream.str());
// }
if(tmodel_opt.empty()){
cout << "Warning: model time sequence is not provided, self generating time sequence from model input" << endl;
@@ -233,21 +233,21 @@ int main(int argc,char **argv) {
}
if(find(direction_opt.begin(),direction_opt.end(),"smooth")!=direction_opt.end()){
if(outputfw_opt.empty()){
- errorStream << "Error: output forward dataset is not provided, use option -ofw" << endl;
- throw(errorStream.str());
+ errorStream << "Error: output forward dataset is not provided, use option -ofw" << endl;
+ throw(errorStream.str());
}
if(outputbw_opt.empty()){
- errorStream << "Error: output backward datasets is not provided, use option -obw" << endl;
- throw(errorStream.str());
+ errorStream << "Error: output backward datasets is not provided, use option -obw" << endl;
+ throw(errorStream.str());
}
if(outputfb_opt.empty()){
- errorStream << "Error: output smooth datasets is not provided, use option -ofb" << endl;
- throw(errorStream.str());
+ errorStream << "Error: output smooth datasets is not provided, use option -ofb" << endl;
+ throw(errorStream.str());
}
if(uncertfb_opt.empty()){
- ostringstream uncertStream;
- uncertStream << outputfb_opt[0] << "_uncert";
- uncertfb_opt.push_back(uncertStream.str());
+ ostringstream uncertStream;
+ uncertStream << outputfb_opt[0] << "_uncert";
+ uncertfb_opt.push_back(uncertStream.str());
}
}
}
@@ -323,7 +323,13 @@ int main(int argc,char **argv) {
if(down_opt.empty()){
double resModel=imgReaderModel1.getDeltaX();
double resObs=imgReaderObs.getDeltaX();
- int down=static_cast<int>(ceil(resModel/resObs));
+ // int down=static_cast<int>(ceil(resModel/resObs));
+ double down_f=resModel/resObs;
+ // round up to the next integer value, unless we already have
+ // an integer (or a value close enough)
+ int down=static_cast<int>(ceil(down_f - eps_opt[0]));
+ // the downsampling factor should be odd
+
if(!(down%2))
down+=1;
down_opt.push_back(down);
@@ -343,13 +349,13 @@ int main(int argc,char **argv) {
try{
if(tobservation_opt.size()<nobs){
if(nobs==nmodel){
- while(tobservation_opt.size()<nobs)
- tobservation_opt.push_back(tobservation_opt.size());
+ while(tobservation_opt.size()<nobs)
+ tobservation_opt.push_back(tobservation_opt.size());
}
else{
- ostringstream errorStream;
- errorStream << "Error: please provide time sequence for observation using option tobs" << endl;
- throw(errorStream.str());
+ ostringstream errorStream;
+ errorStream << "Error: please provide time sequence for observation using option tobs" << endl;
+ throw(errorStream.str());
}
}
}
@@ -357,7 +363,7 @@ int main(int argc,char **argv) {
std::cout << errorString << std::endl;
exit(1);
}
-
+
vector<int> relobsindex;
for(int tindex=0;tindex<tobservation_opt.size();++tindex){
@@ -369,13 +375,13 @@ int main(int argc,char **argv) {
if(verbose_opt[0]){
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();
if(observation_opt.size()>tindex)
- cout << ", filename observation: " << observation_opt[tindex];
+ cout << ", filename observation: " << observation_opt[tindex];
else
- cout << ", observation band index: " << tindex;
+ cout << ", observation band index: " << tindex;
if(model_opt.size()>relpos)
- cout << ", filename of corresponding model: " << model_opt[relpos] << endl;
+ cout << ", filename of corresponding model: " << model_opt[relpos] << endl;
else
- cout << ", band index of corresponding model: " << relpos;
+ cout << ", band index of corresponding model: " << relpos;
}
}
@@ -399,8 +405,8 @@ int main(int argc,char **argv) {
cout << "Running forward model" << endl;
obsindex=0;
if(verbose_opt[0])
- cout << "Opening image " << outputfw_opt[0] << " for writing " << endl << flush;
-
+ cout << "Opening image " << outputfw_opt[0] << " for writing " << endl << flush;
+
ImgWriterGdal imgWriterEst;
ImgWriterGdal imgWriterUncert;
imgWriterEst.open(outputfw_opt[0],ncol,nrow,nmodel,theType,imageType,option_opt);
@@ -412,39 +418,40 @@ int main(int argc,char **argv) {
imgWriterUncert.setGeoTransform(geotransform);
try{
- //test
- if(gain_opt.size()){
- imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
- imgWriterGain.setProjectionProj4(projection_opt[0]);
- imgWriterGain.setGeoTransform(geotransform);
- imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
- }
-
- if(verbose_opt[0]){
- cout << "processing time " << tmodel_opt[0] << endl;
- if(obsindex<relobsindex.size()){
- assert(tmodel_opt.size()>relobsindex[obsindex]);
- cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
- }
- else
- cout << "There is no next observation" << endl;
- }
- if(model_opt.size()==nmodel){
- imgReaderModel1.open(model_opt[0]);
- imgReaderModel1.setNoData(modnodata_opt);
- }
- if(modelmask_opt.size()==nmodel){
- imgReaderModel1Mask.open(modelmask_opt[0]);
- imgReaderModel1Mask.setNoData(msknodata_opt);
- }
+ //test
+ if(gain_opt.size()){
+ imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
+ imgWriterGain.setProjectionProj4(projection_opt[0]);
+ imgWriterGain.setGeoTransform(geotransform);
+ imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
+ }
+
+ if(verbose_opt[0]){
+ cout << "processing time " << tmodel_opt[0] << endl;
+ if(obsindex<relobsindex.size()){
+ assert(tmodel_opt.size()>relobsindex[obsindex]);
+ cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+ }
+ else
+ cout << "There is no next observation" << endl;
+ }
+ if(model_opt.size()==nmodel){
+ imgReaderModel1.open(model_opt[0]);
+ imgReaderModel1.setNoData(modnodata_opt);
+ }
+ if(modelmask_opt.size()==nmodel){
+ imgReaderModel1Mask.open(modelmask_opt[0]);
+ imgReaderModel1Mask.setNoData(msknodata_opt);
+ }
}
catch(string errorString){
- cerr << errorString << endl;
+ cerr << errorString << endl;
}
catch(...){
- cerr << "Error opening file " << model_opt[0] << endl;
+ cerr << "Error opening file " << model_opt[0] << endl;
}
-
+
+ double modEps=0.00001;
double modRow=0;
double modCol=0;
double lowerCol=0;
@@ -452,601 +459,667 @@ int main(int argc,char **argv) {
RESAMPLE theResample=BILINEAR;
if(relobsindex[0]>0){//initialize output_opt[0] as model[0]
- //write first model as output
- if(verbose_opt[0])
- cout << "write first model as output" << endl;
- for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
- vector<double> estReadBuffer;
- vector<double> lineModelMask;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- //test
- vector<double> gainWriteBuffer(ncol);
- try{
- for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
- imgWriterEst.image2geo(0,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
- cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- }
-
- int readModelBand=(model_opt.size()==nmodel)? 0:0;
- int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
- imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
- if(modelmask_opt.size())
- imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
- for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
- for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
- imgWriterEst.image2geo(icol,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- if(modelmask_opt.size()){
- if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- gainWriteBuffer[icol]=obsnodata_opt[0];
- continue;
- }
- }
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
- // double modValue=estReadBuffer[modCol];
- if(imgReaderModel1.isNoData(modValue)){
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- gainWriteBuffer[icol]=obsnodata_opt[0];
- continue;
- }
- estWriteBuffer[icol]=modValue;
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- }
- uncertWriteBuffer[icol]=uncertModel_opt[0];
- //test
- gainWriteBuffer[icol]=0;
- }
- }
- imgWriterEst.writeData(estWriteBuffer,irow,0);
- imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
- //test
- if(gain_opt.size())
- imgWriterGain.writeData(gainWriteBuffer,irow,0);
- }
- }
- catch(string errorString){
- cerr << errorString << endl;
- }
- catch(...){
- cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
- }
- }
+ //write first model as output
+ if(verbose_opt[0])
+ cout << "write first model as output" << endl;
+ for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
+ vector<double> estReadBuffer;
+ vector<double> lineModelMask;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ //test
+ vector<double> gainWriteBuffer(ncol);
+ try{
+ for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
+ imgWriterEst.image2geo(0,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
+ cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ }
+
+ int readModelBand=(model_opt.size()==nmodel)? 0:0;
+ int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
+ imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
+ if(modelmask_opt.size())
+ imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
+ for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
+ for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
+ imgWriterEst.image2geo(icol,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ if(modelmask_opt.size()){
+ if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ gainWriteBuffer[icol]=obsnodata_opt[0];
+ continue;
+ }
+ }
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ double modValue;
+ if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
+ modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ }
+ else{
+ modValue=estReadBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue=estReadBuffer[upperCol];
+ }
+ // double modValue=estReadBuffer[modCol];
+
+ if(imgReaderModel1.isNoData(modValue)){
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ gainWriteBuffer[icol]=obsnodata_opt[0];
+ continue;
+ }
+ estWriteBuffer[icol]=modValue;
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ }
+ uncertWriteBuffer[icol]=uncertModel_opt[0];
+ //test
+ gainWriteBuffer[icol]=0;
+ }
+ }
+ imgWriterEst.writeData(estWriteBuffer,irow,0);
+ imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
+ //test
+ if(gain_opt.size())
+ imgWriterGain.writeData(gainWriteBuffer,irow,0);
+ }
+ }
+ catch(string errorString){
+ cerr << errorString << endl;
+ }
+ catch(...){
+ cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
+ }
+ }
}
else{//we have a measurement
- if(verbose_opt[0])
- cout << "we have a measurement at initial time" << endl;
- if(observation_opt.size()==nobs){
- imgReaderObs.open(observation_opt[0]);
- imgReaderObs.setNoData(obsnodata_opt);
- }
- if(observationmask_opt.size()==nobs){
- imgReaderObsMask.open(observationmask_opt[0]);
- imgReaderObsMask.setNoData(msknodata_opt);
- }
- imgReaderObs.getGeoTransform(geotransform);
-
- vector< vector<double> > obsLineVector(down_opt[0]);
- vector<double> obsLineBuffer;
- vector<double> obsMaskLineBuffer;
- vector<double> modelMaskLineBuffer;
- vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
- vector<double> estReadBuffer;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- vector<double> uncertObsLineBuffer;
- //test
- vector<double> gainWriteBuffer(ncol);
-
- if(verbose_opt[0])
- cout << "initialize obsLineVector" << endl;
- assert(down_opt[0]%2);//window size must be odd
- int readObsBand=(observation_opt.size()==nobs)? 0:0;
- int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:0;
- int readModelBand=(model_opt.size()==nmodel)? 0:0;
- int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
- for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
- if(iline<0)//replicate line 0
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
- else
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
- }
- for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
- for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
- imgWriterEst.image2geo(0,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
- if(modelmask_opt.size())
- imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
- int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
- obsLineVector.erase(obsLineVector.begin());
- imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
- obsLineVector.push_back(obsLineBuffer);
-
- if(observationmask_opt.size())
- imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
-
- for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
- for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
- imgWriterEst.image2geo(icol,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- bool modelIsNoData=false;
- if(modelmask_opt.size())
- modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
- // double modValue=estReadBuffer[modCol];
- double errMod=uncertModel_opt[0];//*stdDev*stdDev;
- modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
- bool obsIsNoData=false;
- if(observationmask_opt.size())
- obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
- obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
- if(modelIsNoData){//model is nodata: retain observation
- if(obsIsNoData){//both model and observation nodata
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- gainWriteBuffer[icol]=obsnodata_opt[0];
- }
- else{
- estWriteBuffer[icol]=obsLineBuffer[icol];
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- }
- uncertWriteBuffer[icol]=uncertObs_opt[0];
- }
- }
- else{//model is valid: calculate estimate from model
- estWriteBuffer[icol]=modValue;
- uncertWriteBuffer[icol]=errMod;//in case observation is not valid
- //test
- gainWriteBuffer[icol]=0;
- }
- //measurement update
- if(!obsIsNoData){
- // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
- double kalmanGain=1;
- int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
- int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
- int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
- int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
- obsWindowBuffer.clear();
- for(int iline=0;iline<obsLineVector.size();++iline){
- for(int isample=minCol;isample<=maxCol;++isample){
- assert(isample<obsLineVector[iline].size());
- obsWindowBuffer.push_back(obsLineVector[iline][isample]);
- }
- }
- if(!modelIsNoData){//model is valid
- statfactory::StatFactory statobs;
- statobs.setNoDataValues(obsnodata_opt);
- double obsMeanValue=0;
- double obsVarValue=0;
- statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
- double difference=0;
- difference=obsMeanValue-modValue;
- // errObs=uncertObs_opt[0]*sqrt(difference*difference);
- errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
- // double errorCovariance=errMod;
- double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
- if(errorCovariance+errObs>eps_opt[0])
- kalmanGain=errorCovariance/(errorCovariance+errObs);
- else
- kalmanGain=1;
- estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- if(kalmanGain>=1)
- kalmanGain=1;
- //test
- gainWriteBuffer[icol]=kalmanGain;
- }
- }
- }
- imgWriterEst.writeData(estWriteBuffer,irow,0);
- imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
- //test
- if(gain_opt.size())
- imgWriterGain.writeData(gainWriteBuffer,irow,0);
- }
- }
- if(observation_opt.size()==nobs)
- imgReaderObs.close();
- if(observationmask_opt.size()==nobs)
- imgReaderObsMask.close();
- ++obsindex;
+ if(verbose_opt[0])
+ cout << "we have a measurement at initial time" << endl;
+ if(observation_opt.size()==nobs){
+ imgReaderObs.open(observation_opt[0]);
+ imgReaderObs.setNoData(obsnodata_opt);
+ }
+ if(observationmask_opt.size()==nobs){
+ imgReaderObsMask.open(observationmask_opt[0]);
+ imgReaderObsMask.setNoData(msknodata_opt);
+ }
+ imgReaderObs.getGeoTransform(geotransform);
+
+ vector< vector<double> > obsLineVector(down_opt[0]);
+ vector<double> obsLineBuffer;
+ vector<double> obsMaskLineBuffer;
+ vector<double> modelMaskLineBuffer;
+ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
+ vector<double> estReadBuffer;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ vector<double> uncertObsLineBuffer;
+ //test
+ vector<double> gainWriteBuffer(ncol);
+
+ if(verbose_opt[0])
+ cout << "initialize obsLineVector" << endl;
+ assert(down_opt[0]%2);//window size must be odd
+ int readObsBand=(observation_opt.size()==nobs)? 0:0;
+ int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:0;
+ int readModelBand=(model_opt.size()==nmodel)? 0:0;
+ int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
+ for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
+ else
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
+ }
+ for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
+ for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
+ imgWriterEst.image2geo(0,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
+ if(modelmask_opt.size())
+ imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+ obsLineVector.erase(obsLineVector.begin());
+ imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
+ obsLineVector.push_back(obsLineBuffer);
+
+ if(observationmask_opt.size())
+ imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
+
+ for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
+ for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
+ imgWriterEst.image2geo(icol,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ bool modelIsNoData=false;
+ if(modelmask_opt.size())
+ modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ double modValue;
+ if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
+ modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ }
+ else{
+ modValue=estReadBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue=estReadBuffer[upperCol];
+ }
+ // double modValue=estReadBuffer[modCol];
+
+ double errMod=uncertModel_opt[0];//*stdDev*stdDev;
+ modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
+ bool obsIsNoData=false;
+ if(observationmask_opt.size())
+ obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
+ obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
+ if(modelIsNoData){//model is nodata: retain observation
+ if(obsIsNoData){//both model and observation nodata
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ gainWriteBuffer[icol]=obsnodata_opt[0];
+ }
+ else{
+ estWriteBuffer[icol]=obsLineBuffer[icol];
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ }
+ uncertWriteBuffer[icol]=uncertObs_opt[0];
+ }
+ }
+ else{//model is valid: calculate estimate from model
+ estWriteBuffer[icol]=modValue;
+ uncertWriteBuffer[icol]=errMod;//in case observation is not valid
+ //test
+ gainWriteBuffer[icol]=0;
+ }
+ //measurement update
+ if(!obsIsNoData){
+ // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
+ double kalmanGain=1;
+ int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
+ int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+ obsWindowBuffer.clear();
+ for(int iline=0;iline<obsLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<obsLineVector[iline].size());
+ obsWindowBuffer.push_back(obsLineVector[iline][isample]);
+ }
+ }
+ if(!modelIsNoData){//model is valid
+ statfactory::StatFactory statobs;
+ statobs.setNoDataValues(obsnodata_opt);
+ double obsMeanValue=0;
+ double obsVarValue=0;
+ statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
+ double difference=0;
+ difference=obsMeanValue-modValue;
+ // errObs=uncertObs_opt[0]*sqrt(difference*difference);
+ errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
+ // double errorCovariance=errMod;
+ double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
+ if(errorCovariance+errObs>eps_opt[0])
+ kalmanGain=errorCovariance/(errorCovariance+errObs);
+ else
+ kalmanGain=1;
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ if(kalmanGain>=1)
+ kalmanGain=1;
+ //test
+ gainWriteBuffer[icol]=kalmanGain;
+ }
+ }
+ }
+ imgWriterEst.writeData(estWriteBuffer,irow,0);
+ imgWriterUncert.writeData(uncertWriteBuffer,irow,0);
+ //test
+ if(gain_opt.size())
+ imgWriterGain.writeData(gainWriteBuffer,irow,0);
+ }
+ }
+ if(observation_opt.size()==nobs)
+ imgReaderObs.close();
+ if(observationmask_opt.size()==nobs)
+ imgReaderObsMask.close();
+ ++obsindex;
}
if(model_opt.size()==nmodel)
- imgReaderModel1.close();
+ imgReaderModel1.close();
if(modelmask_opt.size()==nmodel)
- imgReaderModel1Mask.close();
+ imgReaderModel1Mask.close();
imgWriterEst.close();
imgWriterUncert.close();
ImgUpdaterGdal imgUpdaterEst;
ImgUpdaterGdal imgUpdaterUncert;
for(int modindex=1;modindex<nmodel;++modindex){
- imgUpdaterEst.open(outputfw_opt[0]);
- imgUpdaterEst.setNoData(obsnodata_opt);
- imgUpdaterUncert.open(uncertfw_opt[0]);
- if(verbose_opt[0]){
- cout << "processing time " << tmodel_opt[modindex] << endl;
- if(obsindex<relobsindex.size())
- cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
- else
- cout << "There is no next observation" << endl;
- }
-
- //calculate regression between two subsequence model inputs
- if(model_opt.size()==nmodel){
- imgReaderModel1.open(model_opt[modindex-1]);
- imgReaderModel1.setNoData(modnodata_opt);
- imgReaderModel2.open(model_opt[modindex]);
- imgReaderModel2.setNoData(modnodata_opt);
- }
- if(modelmask_opt.size()==nmodel){
- imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
- imgReaderModel1Mask.setNoData(msknodata_opt);
- imgReaderModel2Mask.open(modelmask_opt[modindex]);
- imgReaderModel2Mask.setNoData(msknodata_opt);
- }
-
- pfnProgress(progress,pszMessage,pProgressArg);
-
- bool update=false;
- if(obsindex<relobsindex.size()){
- update=(relobsindex[obsindex]==modindex);
- }
- if(update){
- if(observation_opt.size()==nobs){
- if(verbose_opt[0])
- cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
- imgReaderObs.open(observation_opt[obsindex]);
- imgReaderObs.getGeoTransform(geotransform);
- imgReaderObs.setNoData(obsnodata_opt);
- }
- if(observationmask_opt.size()==nobs){
- imgReaderObsMask.open(observationmask_opt[obsindex]);
- imgReaderObsMask.setNoData(msknodata_opt);
- }
- }
- //prediction (also to fill cloudy pixels in measurement update mode)
- string input;
- input=outputfw_opt[0];
-
- vector< vector<double> > obsLineVector(down_opt[0]);
- vector<double> obsLineBuffer;
- vector<double> obsMaskLineBuffer;
- vector<double> model1MaskLineBuffer;
- vector<double> model2MaskLineBuffer;
- vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
- vector<double> model1LineBuffer;
- vector<double> model2LineBuffer;
- vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
- vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
- vector<double> uncertObsLineBuffer;
- vector< vector<double> > estLineVector(down_opt[0]);
- vector<double> estLineBuffer;
- vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
- vector<double> uncertReadBuffer;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- //test
- vector<double> gainWriteBuffer(ncol);
-
- int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
- int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
- int readModel1Band=(model_opt.size()==nmodel)? 0:modindex-1;
- int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
- int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex-1;
- int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
-
- //initialize obsLineVector if update
- if(update){
- if(verbose_opt[0])
- cout << "initialize obsLineVector" << endl;
- assert(down_opt[0]%2);//window size must be odd
- for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
- if(iline<0)//replicate line 0
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
- else
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
- }
- }
- //initialize estLineVector
- if(verbose_opt[0])
- cout << "initialize estLineVector" << endl;
- assert(down_opt[0]%2);//window size must be odd
-
- for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
- if(iline<0)//replicate line 0
- imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex-1);
- else
- imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex-1);
- }
- statfactory::StatFactory statobs;
- statobs.setNoDataValues(obsnodata_opt);
- for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
- //todo: read entire window for uncertReadBuffer...
- for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
- imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex-1);
- imgUpdaterUncert.image2geo(0,irow,geox,geoy);
- if(model_opt.size()==nmodel){
- imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
- imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- }
- else{
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
- }
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
-
- if(modelmask_opt.size()){
- imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
- if(modelmask_opt.size()==nmodel)
- imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
- else
- imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
- }
-
- int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
- estLineVector.erase(estLineVector.begin());
- imgUpdaterEst.readData(estLineBuffer,maxRow,modindex-1);
- estLineVector.push_back(estLineBuffer);
- estLineBuffer=estLineVector[down_opt[0]/2];
-
- if(update){
- int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
- obsLineVector.erase(obsLineVector.begin());
- imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
- obsLineVector.push_back(obsLineBuffer);
- obsLineBuffer=obsLineVector[down_opt[0]/2];
-
- if(observationmask_opt.size())
- imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
- }
-
- for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
- for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
- imgUpdaterEst.image2geo(icol,irow,geox,geoy);
- int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
- int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
- int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
- int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
- estWindowBuffer.clear();
- for(int iline=0;iline<estLineVector.size();++iline){
- for(int isample=minCol;isample<=maxCol;++isample){
- assert(isample<estLineVector[iline].size());
- estWindowBuffer.push_back(estLineVector[iline][isample]);
- }
- }
- if(update){
- obsWindowBuffer.clear();
- for(int iline=0;iline<obsLineVector.size();++iline){
- for(int isample=minCol;isample<=maxCol;++isample){
- assert(isample<obsLineVector[iline].size());
- obsWindowBuffer.push_back(obsLineVector[iline][isample]);
- }
- }
- }
- double estValue=estLineBuffer[icol];
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- bool model1IsNoData=false;
- if(modelmask_opt.size())
- model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
- model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
- if(model_opt.size()==nmodel)
- imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
- else
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- bool model2IsNoData=false;
- if(modelmask_opt.size())
- model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
- model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
- bool obsIsNoData=false;
- if(observationmask_opt.size())
- obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
- obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
-
- if(imgUpdaterEst.isNoData(estValue)){
- //we have not found any valid data yet, better here to take the current model value if valid
- if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- gainWriteBuffer[icol]=0;
- }
- else{
- estWriteBuffer[icol]=modValue2;
- uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- //test
- gainWriteBuffer[icol]=0;
- }
- }
- else{//previous estimate is valid
- double estMeanValue=0;
- double estVarValue=0;
- statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
- double nvalid=0;
- //time update
- double processNoiseVariance=processNoise_opt[0]*estVarValue;
- //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
-
- if(model1IsNoData||model2IsNoData){
- estWriteBuffer[icol]=estValue;
- // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
- //todo: check following line if makes sense
- uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
- }
- else{//model is good
- double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
- estWriteBuffer[icol]=estValue*modRatio;
- uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
- }
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- //measurement update
- if(update&&!obsIsNoData){
- double kalmanGain=1;
- if(!model2IsNoData){//model is valid
- statfactory::StatFactory statobs;
- statobs.setNoDataValues(obsnodata_opt);
- double obsMeanValue=0;
- double obsVarValue=0;
- double difference=0;
- statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
- difference=obsMeanValue-modValue2;
- // errObs=uncertObs_opt[0]*sqrt(difference*difference);
- errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
-
- if(errObs<eps_opt[0])
- errObs=eps_opt[0];
- double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
-
- if(errorCovariance+errObs>eps_opt[0])
- kalmanGain=errorCovariance/(errorCovariance+errObs);
- else
- kalmanGain=1;
- estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
- uncertWriteBuffer[icol]*=(1-kalmanGain);
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- if(kalmanGain>=1)
- kalmanGain=1;
- //test
- gainWriteBuffer[icol]=kalmanGain;
- }
- }
- }
-
- //test
- if(gain_opt.size())
- imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
- imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
- imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
- progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
- pfnProgress(progress,pszMessage,pProgressArg);
- }
- }
-
- //must close writers to ensure flush
- imgUpdaterEst.close();
- imgUpdaterUncert.close();
-
- if(update){
- if(observation_opt.size()==nobs)
- imgReaderObs.close();
- if(observationmask_opt.size()==nobs)
- imgReaderObsMask.close();
- ++obsindex;
- }
- if(model_opt.size()==nmodel){
- imgReaderModel1.close();
- imgReaderModel2.close();
- }
- if(modelmask_opt.size()==nmodel){
- imgReaderModel1Mask.close();
- imgReaderModel2Mask.close();
- }
+ imgUpdaterEst.open(outputfw_opt[0]);
+ imgUpdaterEst.setNoData(obsnodata_opt);
+ imgUpdaterUncert.open(uncertfw_opt[0]);
+ if(verbose_opt[0]){
+ cout << "processing time " << tmodel_opt[modindex] << endl;
+ if(obsindex<relobsindex.size())
+ cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+ else
+ cout << "There is no next observation" << endl;
+ }
+
+ //calculate regression between two subsequence model inputs
+ if(model_opt.size()==nmodel){
+ imgReaderModel1.open(model_opt[modindex-1]);
+ imgReaderModel1.setNoData(modnodata_opt);
+ imgReaderModel2.open(model_opt[modindex]);
+ imgReaderModel2.setNoData(modnodata_opt);
+ }
+ if(modelmask_opt.size()==nmodel){
+ imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
+ imgReaderModel1Mask.setNoData(msknodata_opt);
+ imgReaderModel2Mask.open(modelmask_opt[modindex]);
+ imgReaderModel2Mask.setNoData(msknodata_opt);
+ }
+
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ bool update=false;
+ if(obsindex<relobsindex.size()){
+ update=(relobsindex[obsindex]==modindex);
+ }
+ if(update){
+ if(observation_opt.size()==nobs){
+ if(verbose_opt[0])
+ cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
+ imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.getGeoTransform(geotransform);
+ imgReaderObs.setNoData(obsnodata_opt);
+ }
+ if(observationmask_opt.size()==nobs){
+ imgReaderObsMask.open(observationmask_opt[obsindex]);
+ imgReaderObsMask.setNoData(msknodata_opt);
+ }
+ }
+ //prediction (also to fill cloudy pixels in measurement update mode)
+ string input;
+ input=outputfw_opt[0];
+
+ vector< vector<double> > obsLineVector(down_opt[0]);
+ vector<double> obsLineBuffer;
+ vector<double> obsMaskLineBuffer;
+ vector<double> model1MaskLineBuffer;
+ vector<double> model2MaskLineBuffer;
+ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
+ vector<double> model1LineBuffer;
+ vector<double> model2LineBuffer;
+ vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
+ vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
+ vector<double> uncertObsLineBuffer;
+ vector< vector<double> > estLineVector(down_opt[0]);
+ vector<double> estLineBuffer;
+ vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
+ vector<double> uncertReadBuffer;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ //test
+ vector<double> gainWriteBuffer(ncol);
+
+ int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
+ int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
+ int readModel1Band=(model_opt.size()==nmodel)? 0:modindex-1;
+ int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
+ int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex-1;
+ int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
+
+ //initialize obsLineVector if update
+ if(update){
+ if(verbose_opt[0])
+ cout << "initialize obsLineVector" << endl;
+ assert(down_opt[0]%2);//window size must be odd
+ for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
+ else
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
+ }
+ }
+ //initialize estLineVector
+ if(verbose_opt[0])
+ cout << "initialize estLineVector" << endl;
+ assert(down_opt[0]%2);//window size must be odd
+
+ for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex-1);
+ else
+ imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex-1);
+ }
+ statfactory::StatFactory statobs;
+ statobs.setNoDataValues(obsnodata_opt);
+ for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
+ //todo: read entire window for uncertReadBuffer...
+ for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
+ imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex-1);
+ imgUpdaterUncert.image2geo(0,irow,geox,geoy);
+ if(model_opt.size()==nmodel){
+ imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
+ // assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
+ imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ }
+ else{
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
}
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
+
+ if(modelmask_opt.size()){
+ imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
+ if(modelmask_opt.size()==nmodel)
+ imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
+ else
+ imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
+ }
+
+ int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
+ estLineVector.erase(estLineVector.begin());
+ imgUpdaterEst.readData(estLineBuffer,maxRow,modindex-1);
+ estLineVector.push_back(estLineBuffer);
+ estLineBuffer=estLineVector[down_opt[0]/2];
+
+ if(update){
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+ obsLineVector.erase(obsLineVector.begin());
+ imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
+ obsLineVector.push_back(obsLineBuffer);
+ obsLineBuffer=obsLineVector[down_opt[0]/2];
+
+ if(observationmask_opt.size())
+ imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
+ }
+
+ for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
+ for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
+ imgUpdaterEst.image2geo(icol,irow,geox,geoy);
+ int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
+ int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
+ estWindowBuffer.clear();
+ for(int iline=0;iline<estLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<estLineVector[iline].size());
+ estWindowBuffer.push_back(estLineVector[iline][isample]);
+ }
+ }
+ if(update){
+ obsWindowBuffer.clear();
+ for(int iline=0;iline<obsLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<obsLineVector[iline].size());
+ obsWindowBuffer.push_back(obsLineVector[iline][isample]);
+ }
+ }
+ }
+ double estValue=estLineBuffer[icol];
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ bool model1IsNoData=false;
+ if(modelmask_opt.size())
+ model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ double modValue1;
+ if(!imgReaderModel1.isNoData(model1LineBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(model1LineBuffer[upperCol])){
+ modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
+ }
+ else{
+ modValue1=model1LineBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue1=model1LineBuffer[upperCol];
+ }
+ // double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
+ model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
+ if(model_opt.size()==nmodel)
+ imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
+ else
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ bool model2IsNoData=false;
+ if(modelmask_opt.size())
+ model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ // double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
+ double modValue2;
+ if(!imgReaderModel1.isNoData(model2LineBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(model2LineBuffer[upperCol])){
+ modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
+ }
+ else{
+ modValue2=model2LineBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue2=model2LineBuffer[upperCol];
+ }
+ model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
+ bool obsIsNoData=false;
+ if(observationmask_opt.size())
+ obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
+ obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
+
+ if(imgUpdaterEst.isNoData(estValue)){
+ //we have not found any valid data yet, better here to take the current model value if valid
+ if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ gainWriteBuffer[icol]=0;
+ }
+ else{
+ estWriteBuffer[icol]=modValue2;
+ uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ //test
+ gainWriteBuffer[icol]=0;
+ }
+ }
+ else{//previous estimate is valid
+ double estMeanValue=0;
+ double estVarValue=0;
+ statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
+ double nvalid=0;
+ //time update
+ double processNoiseVariance=processNoise_opt[0]*estVarValue;
+ //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
+
+ if(model1IsNoData||model2IsNoData){
+ estWriteBuffer[icol]=estValue;
+ // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
+ //todo: check following line if makes sense
+ uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
+ }
+ else{//model is good
+ double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
+ estWriteBuffer[icol]=estValue*modRatio;
+ uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
+ }
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ //measurement update
+ if(update&&!obsIsNoData){
+ double kalmanGain=1;
+ if(!model2IsNoData){//model is valid
+ statfactory::StatFactory statobs;
+ statobs.setNoDataValues(obsnodata_opt);
+ double obsMeanValue=0;
+ double obsVarValue=0;
+ double difference=0;
+ statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
+ difference=obsMeanValue-modValue2;
+ // errObs=uncertObs_opt[0]*sqrt(difference*difference);
+ errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
+
+ if(errObs<eps_opt[0])
+ errObs=eps_opt[0];
+ double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
+
+ if(errorCovariance+errObs>eps_opt[0])
+ kalmanGain=errorCovariance/(errorCovariance+errObs);
+ else
+ kalmanGain=1;
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
+ uncertWriteBuffer[icol]*=(1-kalmanGain);
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ if(kalmanGain>=1)
+ kalmanGain=1;
+ //test
+ gainWriteBuffer[icol]=kalmanGain;
+ }
+ }
+ }
+
//test
if(gain_opt.size())
- imgWriterGain.close();
+ imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
+ imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
+ imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
+ progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ }
+
+ //must close writers to ensure flush
+ imgUpdaterEst.close();
+ imgUpdaterUncert.close();
+
+ if(update){
+ if(observation_opt.size()==nobs)
+ imgReaderObs.close();
+ if(observationmask_opt.size()==nobs)
+ imgReaderObsMask.close();
+ ++obsindex;
+ }
+ if(model_opt.size()==nmodel){
+ imgReaderModel1.close();
+ imgReaderModel2.close();
+ }
+ if(modelmask_opt.size()==nmodel){
+ imgReaderModel1Mask.close();
+ imgReaderModel2Mask.close();
+ }
+ }
+ //test
+ if(gain_opt.size())
+ imgWriterGain.close();
}
}
catch(string errorString){
@@ -1063,7 +1136,7 @@ int main(int argc,char **argv) {
cout << "Running backward model" << endl;
obsindex=relobsindex.size()-1;
if(verbose_opt[0])
- cout << "Opening image " << outputbw_opt[0] << " for writing " << endl;
+ cout << "Opening image " << outputbw_opt[0] << " for writing " << endl;
ImgWriterGdal imgWriterEst;
ImgWriterGdal imgWriterUncert;
@@ -1076,39 +1149,40 @@ int main(int argc,char **argv) {
imgWriterUncert.setGeoTransform(geotransform);
try{
- // //test
- // if(gain_opt.size()){
- // imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
- // imgWriterGain.setProjectionProj4(projection_opt[0]);
- // imgWriterGain.setGeoTransform(geotransform);
- // imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
- // }
-
- if(verbose_opt[0]){
- cout << "processing time " << tmodel_opt.back() << endl;
- if(obsindex<relobsindex.size()){
- assert(tmodel_opt.size()>relobsindex[obsindex]);
- cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
- }
- else
- cout << "There is no next observation" << endl;
- }
- if(model_opt.size()==nmodel){
- imgReaderModel1.open(model_opt.back());
- imgReaderModel1.setNoData(modnodata_opt);
- }
- if(modelmask_opt.size()==nmodel){
- imgReaderModel1Mask.open(modelmask_opt[0]);
- imgReaderModel1Mask.setNoData(msknodata_opt);
- }
+ // //test
+ // if(gain_opt.size()){
+ // imgWriterGain.open(gain_opt[0],ncol,nrow,nmodel,GDT_Float64,imageType,option_opt);
+ // imgWriterGain.setProjectionProj4(projection_opt[0]);
+ // imgWriterGain.setGeoTransform(geotransform);
+ // imgWriterGain.GDALSetNoDataValue(obsnodata_opt[0]);
+ // }
+
+ if(verbose_opt[0]){
+ cout << "processing time " << tmodel_opt.back() << endl;
+ if(obsindex<relobsindex.size()){
+ assert(tmodel_opt.size()>relobsindex[obsindex]);
+ cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+ }
+ else
+ cout << "There is no next observation" << endl;
+ }
+ if(model_opt.size()==nmodel){
+ imgReaderModel1.open(model_opt.back());
+ imgReaderModel1.setNoData(modnodata_opt);
+ }
+ if(modelmask_opt.size()==nmodel){
+ imgReaderModel1Mask.open(modelmask_opt[0]);
+ imgReaderModel1Mask.setNoData(msknodata_opt);
+ }
}
catch(string errorString){
- cerr << errorString << endl;
+ cerr << errorString << endl;
}
catch(...){
- cerr << "Error opening file " << model_opt[0] << endl;
+ cerr << "Error opening file " << model_opt[0] << endl;
}
+ double modEps=0.00001;
double modRow=0;
double modCol=0;
double lowerCol=0;
@@ -1116,603 +1190,667 @@ int main(int argc,char **argv) {
RESAMPLE theResample=BILINEAR;
if(relobsindex.back()<nmodel-1){//initialize output_opt.back() as last model
- //write last model as output
- if(verbose_opt[0])
- cout << "write last model as output" << endl;
- for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
- vector<double> estReadBuffer;
- vector<double> lineModelMask;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- // //test
- // vector<double> gainWriteBuffer(ncol);
- try{
- for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
- imgWriterEst.image2geo(0,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
- cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- }
- // imgReaderModel1.readData(estReadBuffer,modRow,0,theResample);
- int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
- int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
- imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
- if(modelmask_opt.size())
- imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
- for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
- for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
- imgWriterEst.image2geo(icol,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- if(lineModelMask.size()>modCol){
- if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- // gainWriteBuffer[icol]=obsnodata_opt[0];
- continue;
- }
- }
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
- // double modValue=estReadBuffer[modCol];
- if(imgReaderModel1.isNoData(modValue)){
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- // gainWriteBuffer[icol]=obsnodata_opt[0];
- continue;
- }
- estWriteBuffer[icol]=modValue;
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- }
- uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
- //test
- // gainWriteBuffer[icol]=0;
- }
- }
- imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
- imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
- // //test
- // if(gain_opt.size())
- // imgWriterGain.writeData(gainWriteBuffer,irow,0);
- }
- }
- catch(string errorString){
- cerr << errorString << endl;
- }
- catch(...){
- cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
- }
- }
+ //write last model as output
+ if(verbose_opt[0])
+ cout << "write last model as output" << endl;
+ for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
+ vector<double> estReadBuffer;
+ vector<double> lineModelMask;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ // //test
+ // vector<double> gainWriteBuffer(ncol);
+ try{
+ for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
+ imgWriterEst.image2geo(0,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ if(modRow<0||modRow>=imgReaderModel1.nrOfRow()){
+ cerr << "Error: geo coordinates (" << geox << "," << geoy << ") not covered in model image " << imgReaderModel1.getFileName() << endl;
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ }
+ // imgReaderModel1.readData(estReadBuffer,modRow,0,theResample);
+ int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
+ int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:0;
+ imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
+ if(modelmask_opt.size())
+ imgReaderModel1Mask.readData(lineModelMask,modRow,readModelMaskBand,theResample);
+ for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
+ for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
+ imgWriterEst.image2geo(icol,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ if(lineModelMask.size()>modCol){
+ if(imgReaderModel1Mask.isNoData(lineModelMask[modCol])){
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ // gainWriteBuffer[icol]=obsnodata_opt[0];
+ continue;
+ }
+ }
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ double modValue;
+ if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
+ modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ }
+ else{
+ modValue=estReadBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue=estReadBuffer[upperCol];
+ }
+ // double modValue=estReadBuffer[modCol];
+ if(imgReaderModel1.isNoData(modValue)){
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ // gainWriteBuffer[icol]=obsnodata_opt[0];
+ continue;
+ }
+ estWriteBuffer[icol]=modValue;
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ }
+ uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
+ //test
+ // gainWriteBuffer[icol]=0;
+ }
+ }
+ imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
+ imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
+ // //test
+ // if(gain_opt.size())
+ // imgWriterGain.writeData(gainWriteBuffer,irow,0);
+ }
+ }
+ catch(string errorString){
+ cerr << errorString << endl;
+ }
+ catch(...){
+ cerr << "Error writing file " << imgWriterEst.getFileName() << endl;
+ }
+ }
}
else{//we have a measurement at end time
- if(verbose_opt[0])
- cout << "we have a measurement at end time" << endl;
- if(observation_opt.size()==nobs){
- imgReaderObs.open(observation_opt.back());
- imgReaderObs.setNoData(obsnodata_opt);
- }
- if(observationmask_opt.size()==nobs){
- imgReaderObsMask.open(observationmask_opt.back());
- imgReaderObsMask.setNoData(msknodata_opt);
- }
- imgReaderObs.getGeoTransform(geotransform);
-
- vector< vector<double> > obsLineVector(down_opt[0]);
- vector<double> obsLineBuffer;
- vector<double> obsMaskLineBuffer;
- vector<double> modelMaskLineBuffer;
- vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
- vector<double> estReadBuffer;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- vector<double> uncertObsLineBuffer;
- // //test
- // vector<double> gainWriteBuffer(ncol);
-
- if(verbose_opt[0])
- cout << "initialize obsLineVector" << endl;
- assert(down_opt[0]%2);//window size must be odd
- int readObsBand=(observation_opt.size()==nobs)? 0:nobs-1;
- int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:nobs-1;
- int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
- int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:nmodel-1;
- for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
- if(iline<0)//replicate line 0
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
- else
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
- }
- for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
- for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
- imgWriterEst.image2geo(0,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
- if(modelmask_opt.size())
- imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
- int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
- obsLineVector.erase(obsLineVector.begin());
- imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
- obsLineVector.push_back(obsLineBuffer);
-
- if(observationmask_opt.size())
- imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
-
- for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
- for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
- imgWriterEst.image2geo(icol,irow,geox,geoy);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- bool modelIsNoData=false;
- if(modelmask_opt.size())
- modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
- // double modValue=estReadBuffer[modCol];
- double errMod=uncertModel_opt[0];//*stdDev*stdDev;
- modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
- bool obsIsNoData=false;
- if(observationmask_opt.size())
- obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
- obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
- if(modelIsNoData){//model is nodata: retain observation
- if(obsIsNoData){//both model and observation nodata
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- // gainWriteBuffer[icol]=obsnodata_opt[0];
- }
- else{
- estWriteBuffer[icol]=obsLineBuffer[icol];
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- }
- uncertWriteBuffer[icol]=uncertObs_opt[0];
- }
- }
- else{//model is valid: calculate estimate from model
- estWriteBuffer[icol]=modValue;
- uncertWriteBuffer[icol]=errMod;//in case observation is not valid
- //test
- // gainWriteBuffer[icol]=0;
- }
- //measurement update
- if(!obsIsNoData){
- // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
- double kalmanGain=1;
- int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
- int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
- int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
- int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
- obsWindowBuffer.clear();
- for(int iline=0;iline<obsLineVector.size();++iline){
- for(int isample=minCol;isample<=maxCol;++isample){
- assert(isample<obsLineVector[iline].size());
- obsWindowBuffer.push_back(obsLineVector[iline][isample]);
- }
- }
- if(!modelIsNoData){//model is valid
- statfactory::StatFactory statobs;
- statobs.setNoDataValues(obsnodata_opt);
- double obsMeanValue=0;
- double obsVarValue=0;
- statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
- double difference=0;
- difference=obsMeanValue-modValue;
- // errObs=uncertObs_opt[0]*sqrt(difference*difference);
- errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
- // double errorCovariance=errMod;
- double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
- if(errorCovariance+errObs>eps_opt[0])
- kalmanGain=errorCovariance/(errorCovariance+errObs);
- else
- kalmanGain=1;
- estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- if(kalmanGain>=1)
- kalmanGain=1;
- //test
- // gainWriteBuffer[icol]=kalmanGain;
- }
- }
- }
- imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
- imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
- // //test
- // if(gain_opt.size())
- // imgWriterGain.writeData(gainWriteBuffer,irow,0);
- }
- }
- if(observation_opt.size()==nobs)
- imgReaderObs.close();
- if(observationmask_opt.size()==nobs)
- imgReaderObsMask.close();
- --obsindex;
+ if(verbose_opt[0])
+ cout << "we have a measurement at end time" << endl;
+ if(observation_opt.size()==nobs){
+ imgReaderObs.open(observation_opt.back());
+ imgReaderObs.setNoData(obsnodata_opt);
+ }
+ if(observationmask_opt.size()==nobs){
+ imgReaderObsMask.open(observationmask_opt.back());
+ imgReaderObsMask.setNoData(msknodata_opt);
+ }
+ imgReaderObs.getGeoTransform(geotransform);
+
+ vector< vector<double> > obsLineVector(down_opt[0]);
+ vector<double> obsLineBuffer;
+ vector<double> obsMaskLineBuffer;
+ vector<double> modelMaskLineBuffer;
+ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
+ vector<double> estReadBuffer;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ vector<double> uncertObsLineBuffer;
+ // //test
+ // vector<double> gainWriteBuffer(ncol);
+
+ if(verbose_opt[0])
+ cout << "initialize obsLineVector" << endl;
+ assert(down_opt[0]%2);//window size must be odd
+ int readObsBand=(observation_opt.size()==nobs)? 0:nobs-1;
+ int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:nobs-1;
+ int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
+ int readModelMaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:nmodel-1;
+ for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
+ else
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
+ }
+ for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
+ for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
+ imgWriterEst.image2geo(0,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
+ if(modelmask_opt.size())
+ imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+ obsLineVector.erase(obsLineVector.begin());
+ imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
+ obsLineVector.push_back(obsLineBuffer);
+
+ if(observationmask_opt.size())
+ imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsMaskBand);
+
+ for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
+ for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
+ imgWriterEst.image2geo(icol,irow,geox,geoy);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ bool modelIsNoData=false;
+ if(modelmask_opt.size())
+ modelIsNoData=imgReaderModel1Mask.isNoData(modelMaskLineBuffer[modCol]);
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ // double modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ double modValue;
+ if(!imgReaderModel1.isNoData(estReadBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(estReadBuffer[upperCol])){
+ modValue=(modCol-0.5-lowerCol)*estReadBuffer[upperCol]+(1-modCol+0.5+lowerCol)*estReadBuffer[lowerCol];
+ }
+ else{
+ modValue=estReadBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue=estReadBuffer[upperCol];
+ }
+ // double modValue=estReadBuffer[modCol];
+ double errMod=uncertModel_opt[0];//*stdDev*stdDev;
+ modelIsNoData=modelIsNoData||imgReaderModel1.isNoData(modValue);
+ bool obsIsNoData=false;
+ if(observationmask_opt.size())
+ obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
+ obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
+ if(modelIsNoData){//model is nodata: retain observation
+ if(obsIsNoData){//both model and observation nodata
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ // gainWriteBuffer[icol]=obsnodata_opt[0];
+ }
+ else{
+ estWriteBuffer[icol]=obsLineBuffer[icol];
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ }
+ uncertWriteBuffer[icol]=uncertObs_opt[0];
+ }
+ }
+ else{//model is valid: calculate estimate from model
+ estWriteBuffer[icol]=modValue;
+ uncertWriteBuffer[icol]=errMod;//in case observation is not valid
+ //test
+ // gainWriteBuffer[icol]=0;
+ }
+ //measurement update
+ if(!obsIsNoData){
+ // estWriteBuffer[icol]=estReadBuffer[icol]*modValue1/modValue2
+ double kalmanGain=1;
+ int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
+ int maxCol=(icol+down_opt[0]/2<imgReaderObs.nrOfCol()) ? icol+down_opt[0]/2 : imgReaderObs.nrOfCol()-1;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+ obsWindowBuffer.clear();
+ for(int iline=0;iline<obsLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<obsLineVector[iline].size());
+ obsWindowBuffer.push_back(obsLineVector[iline][isample]);
+ }
+ }
+ if(!modelIsNoData){//model is valid
+ statfactory::StatFactory statobs;
+ statobs.setNoDataValues(obsnodata_opt);
+ double obsMeanValue=0;
+ double obsVarValue=0;
+ statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
+ double difference=0;
+ difference=obsMeanValue-modValue;
+ // errObs=uncertObs_opt[0]*sqrt(difference*difference);
+ errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
+ // double errorCovariance=errMod;
+ double errorCovariance=processNoise_opt[0]*obsVarValue;//assumed initial errorCovariance (P in Kalman equations)
+ if(errorCovariance+errObs>eps_opt[0])
+ kalmanGain=errorCovariance/(errorCovariance+errObs);
+ else
+ kalmanGain=1;
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ if(kalmanGain>=1)
+ kalmanGain=1;
+ //test
+ // gainWriteBuffer[icol]=kalmanGain;
+ }
+ }
+ }
+ imgWriterEst.writeData(estWriteBuffer,irow,nmodel-1);
+ imgWriterUncert.writeData(uncertWriteBuffer,irow,nmodel-1);
+ // //test
+ // if(gain_opt.size())
+ // imgWriterGain.writeData(gainWriteBuffer,irow,0);
+ }
+ }
+ if(observation_opt.size()==nobs)
+ imgReaderObs.close();
+ if(observationmask_opt.size()==nobs)
+ imgReaderObsMask.close();
+ --obsindex;
}
if(model_opt.size()==nmodel)
- imgReaderModel1.close();
+ imgReaderModel1.close();
if(modelmask_opt.size()==nmodel)
- imgReaderModel1Mask.close();
+ imgReaderModel1Mask.close();
imgWriterEst.close();
imgWriterUncert.close();
ImgUpdaterGdal imgUpdaterEst;
ImgUpdaterGdal imgUpdaterUncert;
for(int modindex=nmodel-2;modindex>=0;--modindex){
- imgUpdaterEst.open(outputbw_opt[0]);
- imgUpdaterEst.setNoData(obsnodata_opt);
- imgUpdaterUncert.open(uncertbw_opt[0]);
- if(verbose_opt[0]){
- cout << "processing time " << tmodel_opt[modindex] << endl;
- if(obsindex<relobsindex.size())
- cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
- else
- cout << "There is no next observation" << endl;
- }
-
- //calculate regression between two subsequence model inputs
- if(model_opt.size()==nmodel){
- imgReaderModel1.open(model_opt[modindex+1]);
- imgReaderModel1.setNoData(modnodata_opt);
- imgReaderModel2.open(model_opt[modindex]);
- imgReaderModel2.setNoData(modnodata_opt);
- }
- if(modelmask_opt.size()==nmodel){
- imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
- imgReaderModel1Mask.setNoData(msknodata_opt);
- imgReaderModel2Mask.open(modelmask_opt[modindex]);
- imgReaderModel2Mask.setNoData(msknodata_opt);
- }
-
- pfnProgress(progress,pszMessage,pProgressArg);
-
- bool update=false;
- if(obsindex<relobsindex.size()){
- update=(relobsindex[obsindex]==modindex);
- }
- if(update){
- if(observation_opt.size()==nobs){
- if(verbose_opt[0])
- cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
- imgReaderObs.open(observation_opt[obsindex]);
- imgReaderObs.getGeoTransform(geotransform);
- imgReaderObs.setNoData(obsnodata_opt);
- }
- if(observationmask_opt.size()==nobs){
- imgReaderObsMask.open(observationmask_opt[obsindex]);
- imgReaderObsMask.setNoData(msknodata_opt);
- }
- }
- //prediction (also to fill cloudy pixels in update mode)
- string input;
- input=outputbw_opt[0];
-
- vector< vector<double> > obsLineVector(down_opt[0]);
- vector<double> obsLineBuffer;
- vector<double> obsMaskLineBuffer;
- vector<double> model1MaskLineBuffer;
- vector<double> model2MaskLineBuffer;
- vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
- vector<double> model1LineBuffer;
- vector<double> model2LineBuffer;
- vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
- vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
- vector<double> uncertObsLineBuffer;
- vector< vector<double> > estLineVector(down_opt[0]);
- vector<double> estLineBuffer;
- vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
- vector<double> uncertReadBuffer;
- vector<double> estWriteBuffer(ncol);
- vector<double> uncertWriteBuffer(ncol);
- //test
- // vector<double> gainWriteBuffer(ncol);
-
- int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
- int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
- int readModel1Band=(model_opt.size()==nmodel)? 0:modindex+1;
- int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
- int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex+1;
- int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
-
- //initialize obsLineVector
- if(update){
- if(verbose_opt[0])
- cout << "initialize obsLineVector" << endl;
- assert(down_opt[0]%2);//window size must be odd
- for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
- if(iline<0)//replicate line 0
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
- else
- imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
- }
- }
- //initialize estLineVector
- if(verbose_opt[0])
- cout << "initialize estLineVector" << endl;
- assert(down_opt[0]%2);//window size must be odd
-
- for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
- if(iline<0)//replicate line 0
- imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex+1);
- else
- imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex+1);
- }
- statfactory::StatFactory statobs;
- statobs.setNoDataValues(obsnodata_opt);
-
- for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
- //todo: read entire window for uncertReadBuffer...
- for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
- imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex+1);
- imgUpdaterUncert.image2geo(0,irow,geox,geoy);
- if(model_opt.size()==nmodel){
- imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
- assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
- imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- }
- else{
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
- }
-
- assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
- imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
- if(modelmask_opt.size()){
- imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
- if(modelmask_opt.size()==nmodel)
- imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
- else
- imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
- }
- int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
- estLineVector.erase(estLineVector.begin());
- imgUpdaterEst.readData(estLineBuffer,maxRow,modindex+1);
- estLineVector.push_back(estLineBuffer);
- estLineBuffer=estLineVector[down_opt[0]/2];
-
- if(update){
- int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
- obsLineVector.erase(obsLineVector.begin());
- imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
- obsLineVector.push_back(obsLineBuffer);
- obsLineBuffer=obsLineVector[down_opt[0]/2];
-
- if(observationmask_opt.size())
- imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsBand);
- }
- for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
- for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
- imgUpdaterEst.image2geo(icol,irow,geox,geoy);
- int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
- int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
- int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
- int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
- estWindowBuffer.clear();
- for(int iline=0;iline<estLineVector.size();++iline){
- for(int isample=minCol;isample<=maxCol;++isample){
- assert(isample<estLineVector[iline].size());
- estWindowBuffer.push_back(estLineVector[iline][isample]);
- }
- }
- if(update){
- obsWindowBuffer.clear();
- for(int iline=0;iline<obsLineVector.size();++iline){
- for(int isample=minCol;isample<=maxCol;++isample){
- assert(isample<obsLineVector[iline].size());
- obsWindowBuffer.push_back(obsLineVector[iline][isample]);
- }
- }
- }
-
- double estValue=estLineBuffer[icol];
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- bool model1IsNoData=false;
-
- if(modelmask_opt.size())
- model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
-
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
- model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
- if(model_opt.size()==nmodel)
- imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
- else
- imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
- bool model2IsNoData=false;
-
- if(modelmask_opt.size())
- model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
- lowerCol=modCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=modCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReaderModel1.nrOfCol())
- upperCol=imgReaderModel1.nrOfCol()-1;
- double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
- model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
- bool obsIsNoData=false;
- if(observationmask_opt.size())
- obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
- obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
-
- if(imgUpdaterEst.isNoData(estValue)){
- //we have not found any valid data yet, better here to take the current model value if valid
- if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- //test
- // gainWriteBuffer[icol]=0;
- }
- else{
- estWriteBuffer[icol]=modValue2;
- uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- //test
- // gainWriteBuffer[icol]=0;
- }
- }
- else{//previous estimate is valid
- double estMeanValue=0;
- double estVarValue=0;
- statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
- double nvalid=0;
- //time update
- double processNoiseVariance=processNoise_opt[0]*estVarValue;
- //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
-
- if(model1IsNoData||model2IsNoData){
- estWriteBuffer[icol]=estValue;
- // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
- //todo: check following line if makes sense
- uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
- }
- else{//model is good
- double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
- estWriteBuffer[icol]=estValue*modRatio;
- uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
- }
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- //measurement update
- if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
- double kalmanGain=1;
- if(!imgReaderModel1.isNoData(modValue2)){//model is valid
- statfactory::StatFactory statobs;
- statobs.setNoDataValues(obsnodata_opt);
- double obsMeanValue=0;
- double obsVarValue=0;
- double difference=0;
- statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
- difference=obsMeanValue-modValue2;
- // errObs=uncertObs_opt[0]*sqrt(difference*difference);
- errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
-
- if(errObs<eps_opt[0])
- errObs=eps_opt[0];
- double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
-
- if(errorCovariance+errObs>eps_opt[0])
- kalmanGain=errorCovariance/(errorCovariance+errObs);
- else
- kalmanGain=1;
- estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
- uncertWriteBuffer[icol]*=(1-kalmanGain);
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- if(kalmanGain>=1)
- kalmanGain=1;
- //test
- // gainWriteBuffer[icol]=kalmanGain;
- }
- }
- }
- // //test
- // if(gain_opt.size())
- // imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
- imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
- imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
- progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
- pfnProgress(progress,pszMessage,pProgressArg);
- }
- }
- //must close writers to ensure flush
- imgUpdaterEst.close();
- imgUpdaterUncert.close();
-
- if(update){
- if(observation_opt.size()==nobs)
- imgReaderObs.close();
- if(observationmask_opt.size()==nobs)
- imgReaderObsMask.close();
- --obsindex;
- }
- if(model_opt.size()==nmodel){
- imgReaderModel1.close();
- imgReaderModel2.close();
- }
- if(modelmask_opt.size()==nmodel){
- imgReaderModel1Mask.close();
- imgReaderModel2Mask.close();
- }
+ imgUpdaterEst.open(outputbw_opt[0]);
+ imgUpdaterEst.setNoData(obsnodata_opt);
+ imgUpdaterUncert.open(uncertbw_opt[0]);
+ if(verbose_opt[0]){
+ cout << "processing time " << tmodel_opt[modindex] << endl;
+ if(obsindex<relobsindex.size())
+ cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+ else
+ cout << "There is no next observation" << endl;
+ }
+
+ //calculate regression between two subsequence model inputs
+ if(model_opt.size()==nmodel){
+ imgReaderModel1.open(model_opt[modindex+1]);
+ imgReaderModel1.setNoData(modnodata_opt);
+ imgReaderModel2.open(model_opt[modindex]);
+ imgReaderModel2.setNoData(modnodata_opt);
+ }
+ if(modelmask_opt.size()==nmodel){
+ imgReaderModel1Mask.open(modelmask_opt[modindex-1]);
+ imgReaderModel1Mask.setNoData(msknodata_opt);
+ imgReaderModel2Mask.open(modelmask_opt[modindex]);
+ imgReaderModel2Mask.setNoData(msknodata_opt);
+ }
+
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ bool update=false;
+ if(obsindex<relobsindex.size()){
+ update=(relobsindex[obsindex]==modindex);
+ }
+ if(update){
+ if(observation_opt.size()==nobs){
+ if(verbose_opt[0])
+ cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
+ imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.getGeoTransform(geotransform);
+ imgReaderObs.setNoData(obsnodata_opt);
+ }
+ if(observationmask_opt.size()==nobs){
+ imgReaderObsMask.open(observationmask_opt[obsindex]);
+ imgReaderObsMask.setNoData(msknodata_opt);
+ }
+ }
+ //prediction (also to fill cloudy pixels in update mode)
+ string input;
+ input=outputbw_opt[0];
+
+ vector< vector<double> > obsLineVector(down_opt[0]);
+ vector<double> obsLineBuffer;
+ vector<double> obsMaskLineBuffer;
+ vector<double> model1MaskLineBuffer;
+ vector<double> model2MaskLineBuffer;
+ vector<double> obsWindowBuffer;//buffer for observation to calculate average corresponding to model pixel
+ vector<double> model1LineBuffer;
+ vector<double> model2LineBuffer;
+ vector<double> model1buffer;//buffer for model 1 to calculate time regression based on window
+ vector<double> model2buffer;//buffer for model 2 to calculate time regression based on window
+ vector<double> uncertObsLineBuffer;
+ vector< vector<double> > estLineVector(down_opt[0]);
+ vector<double> estLineBuffer;
+ vector<double> estWindowBuffer;//buffer for estimate to calculate average corresponding to model pixel
+ vector<double> uncertReadBuffer;
+ vector<double> estWriteBuffer(ncol);
+ vector<double> uncertWriteBuffer(ncol);
+ //test
+ // vector<double> gainWriteBuffer(ncol);
+
+ int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
+ int readObsMaskBand=(observationmask_opt.size()==nobs)? mskband_opt[0]:obsindex;
+ int readModel1Band=(model_opt.size()==nmodel)? 0:modindex+1;
+ int readModel2Band=(model_opt.size()==nmodel)? 0:modindex;
+ int readModel1MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex+1;
+ int readModel2MaskBand=(modelmask_opt.size()==nmodel)? mskband_opt[0]:modindex;
+
+ //initialize obsLineVector
+ if(update){
+ if(verbose_opt[0])
+ cout << "initialize obsLineVector" << endl;
+ assert(down_opt[0]%2);//window size must be odd
+ for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],0,readObsBand);
+ else
+ imgReaderObs.readData(obsLineVector[iline+down_opt[0]/2],iline,readObsBand);
+ }
+ }
+ //initialize estLineVector
+ if(verbose_opt[0])
+ cout << "initialize estLineVector" << endl;
+ assert(down_opt[0]%2);//window size must be odd
+
+ for(int iline=-down_opt[0]/2;iline<down_opt[0]/2+1;++iline){
+ if(iline<0)//replicate line 0
+ imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],0,modindex+1);
+ else
+ imgUpdaterEst.readData(estLineVector[iline+down_opt[0]/2],iline,modindex+1);
+ }
+ statfactory::StatFactory statobs;
+ statobs.setNoDataValues(obsnodata_opt);
+
+ for(int jrow=0;jrow<nrow;jrow+=down_opt[0]){
+ //todo: read entire window for uncertReadBuffer...
+ for(int irow=jrow;irow<jrow+down_opt[0]&&irow<nrow;++irow){
+ imgUpdaterUncert.readData(uncertReadBuffer,irow,modindex+1);
+ imgUpdaterUncert.image2geo(0,irow,geox,geoy);
+ if(model_opt.size()==nmodel){
+ imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
+ // assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
+ imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ }
+ else{
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
+ }
+
+ // assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
+ imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
+ if(modelmask_opt.size()){
+ imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
+ if(modelmask_opt.size()==nmodel)
+ imgReaderModel2Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
+ else
+ imgReaderModel1Mask.readData(model2MaskLineBuffer,modRow,readModel2MaskBand);
+ }
+ int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
+ estLineVector.erase(estLineVector.begin());
+ imgUpdaterEst.readData(estLineBuffer,maxRow,modindex+1);
+ estLineVector.push_back(estLineBuffer);
+ estLineBuffer=estLineVector[down_opt[0]/2];
+
+ if(update){
+ int maxRow=(irow+down_opt[0]/2<imgReaderObs.nrOfRow()) ? irow+down_opt[0]/2 : imgReaderObs.nrOfRow()-1;
+ obsLineVector.erase(obsLineVector.begin());
+ imgReaderObs.readData(obsLineBuffer,maxRow,readObsBand);
+ obsLineVector.push_back(obsLineBuffer);
+ obsLineBuffer=obsLineVector[down_opt[0]/2];
+
+ if(observationmask_opt.size())
+ imgReaderObsMask.readData(obsMaskLineBuffer,irow,readObsBand);
+ }
+ for(int jcol=0;jcol<ncol;jcol+=down_opt[0]){
+ for(int icol=jcol;icol<jcol+down_opt[0]&&icol<ncol;++icol){
+ imgUpdaterEst.image2geo(icol,irow,geox,geoy);
+ int minCol=(icol>down_opt[0]/2) ? icol-down_opt[0]/2 : 0;
+ int maxCol=(icol+down_opt[0]/2<imgUpdaterEst.nrOfCol()) ? icol+down_opt[0]/2 : imgUpdaterEst.nrOfCol()-1;
+ int minRow=(irow>down_opt[0]/2) ? irow-down_opt[0]/2 : 0;
+ int maxRow=(irow+down_opt[0]/2<imgUpdaterEst.nrOfRow()) ? irow+down_opt[0]/2 : imgUpdaterEst.nrOfRow()-1;
+ estWindowBuffer.clear();
+ for(int iline=0;iline<estLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<estLineVector[iline].size());
+ estWindowBuffer.push_back(estLineVector[iline][isample]);
+ }
+ }
+ if(update){
+ obsWindowBuffer.clear();
+ for(int iline=0;iline<obsLineVector.size();++iline){
+ for(int isample=minCol;isample<=maxCol;++isample){
+ assert(isample<obsLineVector[iline].size());
+ obsWindowBuffer.push_back(obsLineVector[iline][isample]);
+ }
+ }
+ }
+
+ double estValue=estLineBuffer[icol];
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ bool model1IsNoData=false;
+
+ if(modelmask_opt.size())
+ model1IsNoData=imgReaderModel1Mask.isNoData(model1MaskLineBuffer[modCol]);
+
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ // double modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
+ double modValue1;
+ if(!imgReaderModel1.isNoData(model1LineBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(model1LineBuffer[upperCol])){
+ modValue1=(modCol-0.5-lowerCol)*model1LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model1LineBuffer[lowerCol];
+ }
+ else{
+ modValue1=model1LineBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue1=model1LineBuffer[upperCol];
+ }
+ model1IsNoData=model1IsNoData||imgReaderModel1.isNoData(modValue1);
+ if(model_opt.size()==nmodel)
+ imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
+ else
+ imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
+ bool model2IsNoData=false;
+
+ if(modelmask_opt.size())
+ model2IsNoData=imgReaderModel1Mask.isNoData(model2MaskLineBuffer[modCol]);
+ // lowerCol=modCol-0.5;
+ // lowerCol=static_cast<int>(lowerCol);
+ // upperCol=modCol+0.5;
+ // upperCol=static_cast<int>(upperCol);
+ lowerCol=static_cast<int>(modCol-0.5+modEps);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(lowerCol>=imgReaderModel1.nrOfCol())
+ lowerCol=imgReaderModel1.nrOfCol()-1;
+ upperCol=lowerCol+1.0;
+ if(upperCol>=imgReaderModel1.nrOfCol())
+ upperCol=imgReaderModel1.nrOfCol()-1;
+ double modValue2;
+ if(!imgReaderModel1.isNoData(model2LineBuffer[lowerCol])){
+ if(!imgReaderModel1.isNoData(model2LineBuffer[upperCol])){
+ modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
+ }
+ else{
+ modValue2=model2LineBuffer[lowerCol];
+ }
+ }
+ else{
+ modValue2=model2LineBuffer[upperCol];
+ }
+ // double modValue2=(modCol-0.5-lowerCol)*model2LineBuffer[upperCol]+(1-modCol+0.5+lowerCol)*model2LineBuffer[lowerCol];
+ model2IsNoData=model2IsNoData||imgReaderModel1.isNoData(modValue2);
+ bool obsIsNoData=false;
+ if(observationmask_opt.size())
+ obsIsNoData=imgReaderObsMask.isNoData(obsMaskLineBuffer[icol]);
+ obsIsNoData=obsIsNoData||imgReaderObs.isNoData(obsLineBuffer[icol]);
+
+ if(imgUpdaterEst.isNoData(estValue)){
+ //we have not found any valid data yet, better here to take the current model value if valid
+ if(model2IsNoData){//if both estimate and model are no-data, set obs to nodata
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ //test
+ // gainWriteBuffer[icol]=0;
+ }
+ else{
+ estWriteBuffer[icol]=modValue2;
+ uncertWriteBuffer[icol]=uncertModel_opt[0];//*stdDev*stdDev;
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ //test
+ // gainWriteBuffer[icol]=0;
+ }
+ }
+ else{//previous estimate is valid
+ double estMeanValue=0;
+ double estVarValue=0;
+ statobs.meanVar(estWindowBuffer,estMeanValue,estVarValue);
+ double nvalid=0;
+ //time update
+ double processNoiseVariance=processNoise_opt[0]*estVarValue;
+ //estimate stability of weight distribution from model (low resolution) data in a window mod1 -> mod2 and assume distribution holds at fine spatial resolution.
+
+ if(model1IsNoData||model2IsNoData){
+ estWriteBuffer[icol]=estValue;
+ // uncertWriteBuffer[icol]=uncertReadBuffer[icol]+processNoiseVariance;
+ //todo: check following line if makes sense
+ uncertWriteBuffer[icol]=uncertReadBuffer[icol]+uncertObs_opt[0];
+ }
+ else{//model is good
+ double modRatio=modValue2/modValue1;//transition matrix A in Kalman equations
+ estWriteBuffer[icol]=estValue*modRatio;
+ uncertWriteBuffer[icol]=uncertReadBuffer[icol]*modRatio*modRatio+processNoiseVariance;
+ }
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ //measurement update
+ if(update&&!imgReaderObs.isNoData(obsLineBuffer[icol])){
+ double kalmanGain=1;
+ if(!imgReaderModel1.isNoData(modValue2)){//model is valid
+ statfactory::StatFactory statobs;
+ statobs.setNoDataValues(obsnodata_opt);
+ double obsMeanValue=0;
+ double obsVarValue=0;
+ double difference=0;
+ statobs.meanVar(obsWindowBuffer,obsMeanValue,obsVarValue);
+ difference=obsMeanValue-modValue2;
+ // errObs=uncertObs_opt[0]*sqrt(difference*difference);
+ errObs=uncertObs_opt[0]*difference*difference;//uncertainty of the observation (R in Kalman equations)
+
+ if(errObs<eps_opt[0])
+ errObs=eps_opt[0];
+ double errorCovariance=uncertWriteBuffer[icol];//P in Kalman equations
+
+ if(errorCovariance+errObs>eps_opt[0])
+ kalmanGain=errorCovariance/(errorCovariance+errObs);
+ else
+ kalmanGain=1;
+ estWriteBuffer[icol]+=kalmanGain*(obsLineBuffer[icol]-estWriteBuffer[icol]);
+ uncertWriteBuffer[icol]*=(1-kalmanGain);
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ if(kalmanGain>=1)
+ kalmanGain=1;
+ //test
+ // gainWriteBuffer[icol]=kalmanGain;
+ }
+ }
}
// //test
// if(gain_opt.size())
- // imgWriterGain.close();
+ // imgWriterGain.writeData(gainWriteBuffer,irow,modindex);
+ imgUpdaterEst.writeData(estWriteBuffer,irow,modindex);
+ imgUpdaterUncert.writeData(uncertWriteBuffer,irow,modindex);
+ progress=static_cast<float>((irow+1.0)/imgUpdaterEst.nrOfRow());
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ }
+ //must close writers to ensure flush
+ imgUpdaterEst.close();
+ imgUpdaterUncert.close();
+
+ if(update){
+ if(observation_opt.size()==nobs)
+ imgReaderObs.close();
+ if(observationmask_opt.size()==nobs)
+ imgReaderObsMask.close();
+ --obsindex;
+ }
+ if(model_opt.size()==nmodel){
+ imgReaderModel1.close();
+ imgReaderModel2.close();
+ }
+ if(modelmask_opt.size()==nmodel){
+ imgReaderModel1Mask.close();
+ imgReaderModel2Mask.close();
+ }
+ }
+ // //test
+ // if(gain_opt.size())
+ // imgWriterGain.close();
}
}
catch(string errorString){
@@ -1734,7 +1872,7 @@ int main(int argc,char **argv) {
ImgReaderGdal imgReaderBackwardUncert(uncertbw_opt[0]);
imgReaderForward.setNoData(obsnodata_opt);
imgReaderBackward.setNoData(obsnodata_opt);
-
+
assert(imgReaderForward.nrOfBand()==nmodel);
assert(imgReaderForwardUncert.nrOfBand()==nmodel);
assert(imgReaderBackward.nrOfBand()==nmodel);
@@ -1752,11 +1890,11 @@ int main(int argc,char **argv) {
imgWriterUncert.setGeoTransform(geotransform);
for(int modindex=0;modindex<nmodel;++modindex){
if(verbose_opt[0]){
- cout << "processing time " << tmodel_opt[modindex] << endl;
- if(obsindex<relobsindex.size())
- cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
- else
- cout << "There is no next observation" << endl;
+ cout << "processing time " << tmodel_opt[modindex] << endl;
+ if(obsindex<relobsindex.size())
+ cout << "next observation " << tmodel_opt[relobsindex[obsindex]] << endl;
+ else
+ cout << "There is no next observation" << endl;
}
vector<double> estForwardBuffer;
vector<double> estBackwardBuffer;
@@ -1769,94 +1907,94 @@ int main(int argc,char **argv) {
bool update=false;
if(obsindex<relobsindex.size()){
- update=(relobsindex[obsindex]==modindex);
+ update=(relobsindex[obsindex]==modindex);
}
int readObsBand=(observation_opt.size()==nobs)? 0:obsindex;
if(update){
- if(observation_opt.size()==nobs){
- if(verbose_opt[0])
- cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
- imgReaderObs.open(observation_opt[obsindex]);
- imgReaderObs.setNoData(obsnodata_opt);
- imgReaderObs.getGeoTransform(geotransform);
- }
- if(observationmask_opt.size()==nobs){
- imgReaderObsMask.open(observationmask_opt[obsindex]);
- imgReaderObsMask.setNoData(msknodata_opt);
- }
+ if(observation_opt.size()==nobs){
+ if(verbose_opt[0])
+ cout << "***update " << relobsindex[obsindex] << " = " << modindex << " " << observation_opt[obsindex] << " ***" << endl;
+ imgReaderObs.open(observation_opt[obsindex]);
+ imgReaderObs.setNoData(obsnodata_opt);
+ imgReaderObs.getGeoTransform(geotransform);
+ }
+ if(observationmask_opt.size()==nobs){
+ imgReaderObsMask.open(observationmask_opt[obsindex]);
+ imgReaderObsMask.setNoData(msknodata_opt);
+ }
}
pfnProgress(progress,pszMessage,pProgressArg);
for(int irow=0;irow<imgWriterEst.nrOfRow();++irow){
- assert(irow<imgReaderForward.nrOfRow());
- assert(irow<imgReaderBackward.nrOfRow());
- imgReaderForward.readData(estForwardBuffer,irow,modindex);
- imgReaderBackward.readData(estBackwardBuffer,irow,modindex);
- imgReaderForwardUncert.readData(uncertForwardBuffer,irow,modindex);
- imgReaderBackwardUncert.readData(uncertBackwardBuffer,irow,modindex);
-
- if(update){
- if(observation_opt.size()==nobs)
- imgReaderObs.readData(estWriteBuffer,irow,readObsBand);
- if(observationmask_opt.size())
- imgReaderObsMask.readData(uncertObsLineBuffer,irow,readObsBand);
- }
-
- for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
- imgWriterEst.image2geo(icol,irow,geox,geoy);
- double A=estForwardBuffer[icol];
- double B=estBackwardBuffer[icol];
- double C=uncertForwardBuffer[icol];
- double D=uncertBackwardBuffer[icol];
- double uncertObs=uncertObs_opt[0];
-
- double noemer=(C+D);
- //todo: consistently check for division by zero...
- if(imgReaderForward.isNoData(A)&&imgReaderBackward.isNoData(B)){
- estWriteBuffer[icol]=obsnodata_opt[0];
- uncertWriteBuffer[icol]=uncertNodata_opt[0];
- }
- else if(imgReaderForward.isNoData(A)){
- estWriteBuffer[icol]=B;
- uncertWriteBuffer[icol]=uncertBackwardBuffer[icol];
- }
- else if(imgReaderForward.isNoData(B)){
- estWriteBuffer[icol]=A;
- uncertWriteBuffer[icol]=uncertForwardBuffer[icol];
- }
- else{
- if(noemer<eps_opt[0]){//simple average if both uncertainties are ~>0
- estWriteBuffer[icol]=0.5*(A+B);
- uncertWriteBuffer[icol]=eps_opt[0];
- }
- else{
- estWriteBuffer[icol]=(A*D+B*C)/noemer;
- uncertWriteBuffer[icol]=C*D/noemer;
- if(obsmin_opt.size()){
- if(estWriteBuffer[icol]<obsmin_opt[0])
- estWriteBuffer[icol]=obsmin_opt[0];
- }
- if(obsmax_opt.size()){
- if(estWriteBuffer[icol]>obsmax_opt[0])
- estWriteBuffer[icol]=obsmax_opt[0];
- if(uncertWriteBuffer[icol]>obsmax_opt[0])
- uncertWriteBuffer[icol]=obsmax_opt[0];
- }
- }
- }
- }
- imgWriterEst.writeData(estWriteBuffer,irow,modindex);
- imgWriterUncert.writeData(uncertWriteBuffer,irow,modindex);
- progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
- pfnProgress(progress,pszMessage,pProgressArg);
+ assert(irow<imgReaderForward.nrOfRow());
+ assert(irow<imgReaderBackward.nrOfRow());
+ imgReaderForward.readData(estForwardBuffer,irow,modindex);
+ imgReaderBackward.readData(estBackwardBuffer,irow,modindex);
+ imgReaderForwardUncert.readData(uncertForwardBuffer,irow,modindex);
+ imgReaderBackwardUncert.readData(uncertBackwardBuffer,irow,modindex);
+
+ if(update){
+ if(observation_opt.size()==nobs)
+ imgReaderObs.readData(estWriteBuffer,irow,readObsBand);
+ if(observationmask_opt.size())
+ imgReaderObsMask.readData(uncertObsLineBuffer,irow,readObsBand);
+ }
+
+ for(int icol=0;icol<imgWriterEst.nrOfCol();++icol){
+ imgWriterEst.image2geo(icol,irow,geox,geoy);
+ double A=estForwardBuffer[icol];
+ double B=estBackwardBuffer[icol];
+ double C=uncertForwardBuffer[icol];
+ double D=uncertBackwardBuffer[icol];
+ double uncertObs=uncertObs_opt[0];
+
+ double noemer=(C+D);
+ //todo: consistently check for division by zero...
+ if(imgReaderForward.isNoData(A)&&imgReaderBackward.isNoData(B)){
+ estWriteBuffer[icol]=obsnodata_opt[0];
+ uncertWriteBuffer[icol]=uncertNodata_opt[0];
+ }
+ else if(imgReaderForward.isNoData(A)){
+ estWriteBuffer[icol]=B;
+ uncertWriteBuffer[icol]=uncertBackwardBuffer[icol];
+ }
+ else if(imgReaderForward.isNoData(B)){
+ estWriteBuffer[icol]=A;
+ uncertWriteBuffer[icol]=uncertForwardBuffer[icol];
+ }
+ else{
+ if(noemer<eps_opt[0]){//simple average if both uncertainties are ~>0
+ estWriteBuffer[icol]=0.5*(A+B);
+ uncertWriteBuffer[icol]=eps_opt[0];
+ }
+ else{
+ estWriteBuffer[icol]=(A*D+B*C)/noemer;
+ uncertWriteBuffer[icol]=C*D/noemer;
+ if(obsmin_opt.size()){
+ if(estWriteBuffer[icol]<obsmin_opt[0])
+ estWriteBuffer[icol]=obsmin_opt[0];
+ }
+ if(obsmax_opt.size()){
+ if(estWriteBuffer[icol]>obsmax_opt[0])
+ estWriteBuffer[icol]=obsmax_opt[0];
+ if(uncertWriteBuffer[icol]>obsmax_opt[0])
+ uncertWriteBuffer[icol]=obsmax_opt[0];
+ }
+ }
+ }
+ }
+ imgWriterEst.writeData(estWriteBuffer,irow,modindex);
+ imgWriterUncert.writeData(uncertWriteBuffer,irow,modindex);
+ progress=static_cast<float>((irow+1.0)/imgWriterEst.nrOfRow());
+ pfnProgress(progress,pszMessage,pProgressArg);
}
if(update){
- if(observation_opt.size()==nobs)
- imgReaderObs.close();
- ++obsindex;
+ if(observation_opt.size()==nobs)
+ imgReaderObs.close();
+ ++obsindex;
}
}
imgReaderForward.close();
@@ -1869,4 +2007,3 @@ int main(int argc,char **argv) {
if(model_opt.size()<nmodel)
imgReaderModel1.close();
}
-
diff --git a/src/apps/pkstat.cc b/src/apps/pkstat.cc
index c2af39f..b661316 100644
--- a/src/apps/pkstat.cc
+++ b/src/apps/pkstat.cc
@@ -25,52 +25,52 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "algorithms/ImgRegression.h"
/******************************************************************************/
/*! \page pkstat pkstat
- program to calculate basic statistics from raster dataset
-## SYNOPSIS
+ program to calculate basic statistics from raster dataset
+ ## SYNOPSIS
+
+ <code>
+ Usage: pkstat -i input
+ </code>
+
+ \section pkstat_options Options
+ - use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
+ - short option `-h` shows basic options only, long option `--help` shows all options
+ |short|long|type|default|description|
+ |-----|----|----|-------|-----------|
+ | i | input | std::string | |name of the input raster dataset |
+ | b | band | unsigned short | 0 |band(s) on which to calculate statistics |
+ | f | filename | bool | false |Shows image filename |
+ | stats | statistics | bool | false |Shows basic statistics (min,max, mean and stdDev of the raster datasets) |
+ | nodata | nodata | double | |Set nodata value(s) |
+ | mean | mean | bool | false |calculate mean |
+ | median | median | bool | false |calculate median |
+ | var | var | bool | false |calculate variance |
+ | stdev | stdev | bool | false |calculate standard deviation |
+ | mm | minmax | bool | false |calculate minimum and maximum value |
+ | min | min | bool | false |calculate minimum value |
+ | max | max | bool | false |calculate maximum value |
+ | hist | hist | bool | false |calculate histogram |
+ | nbin | nbin | short | |number of bins to calculate histogram |
+ | rel | relative | bool | false |use percentiles for histogram to calculate histogram |
+ | hist2d | hist2d | bool | false |calculate 2-dimensional histogram based on two images |
+ | cor | correlation | bool | false |calculate Pearson produc-moment correlation coefficient between two raster datasets (defined by -c <col1> -c <col2>) |
+ | rmse | rmse | bool | false |calculate root mean square error between two raster datasets |
+ | reg | regression | bool | false |calculate linear regression between two raster datasets and get correlation coefficient |
+ | regerr | regerr | bool | false |calculate linear regression between two raster datasets and get root mean square error |
+ | preg | preg | bool | false |calculate perpendicular regression between two raster datasets and get correlation coefficient |
+ | ulx | ulx | double | |Upper left x value bounding box |
+ | uly | uly | double | |Upper left y value bounding box |
+ | lrx | lrx | double | |Lower right x value bounding box |
+ | lry | lry | double | |Lower right y value bounding box |
+ | down | down | short | 1 |Down sampling factor (for raster sample datasets only). Can be used to create grid points |
+ | rnd | rnd | unsigned int | 0 |generate random numbers |
+ | scale | scale | double | |Scale(s) for reading input image(s) |
+ | offset | offset | double | |Offset(s) for reading input image(s) |
+ | src_min | src_min | double | |start reading source from this minimum value |
+ | src_max | src_max | double | |stop reading source from this maximum value |
+ | kde | kde | bool | false |Use Kernel density estimation when producing histogram. The standard deviation is estimated based on Silverman's rule of thumb |
-<code>
Usage: pkstat -i input
-</code>
-
-\section pkstat_options Options
- - use either `-short` or `--long` options (both `--long=value` and `--long value` are supported)
- - short option `-h` shows basic options only, long option `--help` shows all options
-|short|long|type|default|description|
-|-----|----|----|-------|-----------|
- | i | input | std::string | |name of the input raster dataset |
- | b | band | unsigned short | 0 |band(s) on which to calculate statistics |
- | f | filename | bool | false |Shows image filename |
- | stats | statistics | bool | false |Shows basic statistics (min,max, mean and stdDev of the raster datasets) |
- | nodata | nodata | double | |Set nodata value(s) |
- | mean | mean | bool | false |calculate mean |
- | median | median | bool | false |calculate median |
- | var | var | bool | false |calculate variance |
- | stdev | stdev | bool | false |calculate standard deviation |
- | mm | minmax | bool | false |calculate minimum and maximum value |
- | min | min | bool | false |calculate minimum value |
- | max | max | bool | false |calculate maximum value |
- | hist | hist | bool | false |calculate histogram |
- | nbin | nbin | short | |number of bins to calculate histogram |
- | rel | relative | bool | false |use percentiles for histogram to calculate histogram |
- | hist2d | hist2d | bool | false |calculate 2-dimensional histogram based on two images |
- | cor | correlation | bool | false |calculate Pearson produc-moment correlation coefficient between two raster datasets (defined by -c <col1> -c <col2>) |
- | rmse | rmse | bool | false |calculate root mean square error between two raster datasets |
- | reg | regression | bool | false |calculate linear regression between two raster datasets and get correlation coefficient |
- | regerr | regerr | bool | false |calculate linear regression between two raster datasets and get root mean square error |
- | preg | preg | bool | false |calculate perpendicular regression between two raster datasets and get correlation coefficient |
- | ulx | ulx | double | |Upper left x value bounding box |
- | uly | uly | double | |Upper left y value bounding box |
- | lrx | lrx | double | |Lower right x value bounding box |
- | lry | lry | double | |Lower right y value bounding box |
- | down | down | short | 1 |Down sampling factor (for raster sample datasets only). Can be used to create grid points |
- | rnd | rnd | unsigned int | 0 |generate random numbers |
- | scale | scale | double | |Scale(s) for reading input image(s) |
- | offset | offset | double | |Offset(s) for reading input image(s) |
- | src_min | src_min | double | |start reading source from this minimum value |
- | src_max | src_max | double | |stop reading source from this maximum value |
- | kde | kde | bool | false |Use Kernel density estimation when producing histogram. The standard deviation is estimated based on Silverman's rule of thumb |
-
-Usage: pkstat -i input
**/
@@ -230,7 +230,7 @@ int main(int argc, char *argv[])
}
if(input_opt.empty()){
std::cerr << "No image dataset provided (use option -i). Use --help for help information";
- exit(0);
+ exit(0);
}
for(int ifile=0;ifile<input_opt.size();++ifile){
try{
@@ -252,8 +252,8 @@ int main(int argc, char *argv[])
minValue=(src_min_opt.size())? src_min_opt[0] : 0;
maxValue=(src_max_opt.size())? src_max_opt[0] : 0;
for(int inodata=0;inodata<nodata_opt.size();++inodata){
- if(!inodata)
- imgReader.GDALSetNoDataValue(nodata_opt[0],band_opt[iband]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
+ if(!inodata)
+ imgReader.GDALSetNoDataValue(nodata_opt[0],band_opt[iband]);//only single no data can be set in GDALRasterBand (used for ComputeStatistics)
}
if(offset_opt.size()>ifile)
@@ -262,91 +262,91 @@ int main(int argc, char *argv[])
imgReader.setScale(scale_opt[ifile],band_opt[iband]);
if(stat_opt[0]||mean_opt[0]||median_opt[0]||var_opt[0]||stdev_opt[0]){//the hard way (in memory)
- statfactory::StatFactory stat;
- vector<double> readBuffer;
- double varValue;
- imgReader.readDataBlock(readBuffer, 0, imgReader.nrOfCol()-1, 0, imgReader.nrOfRow()-1, band_opt[0]);
- stat.setNoDataValues(nodata_opt);
- stat.meanVar(readBuffer,meanValue,varValue);
- medianValue=stat.median(readBuffer);
- stat.minmax(readBuffer,readBuffer.begin(),readBuffer.end(),minValue,maxValue);
- if(mean_opt[0])
- std::cout << "--mean " << meanValue << " ";
- if(median_opt[0])
- std::cout << "--median " << medianValue << " ";
- if(stdev_opt[0])
- std::cout << "--stdDev " << sqrt(varValue) << " ";
- if(var_opt[0])
- std::cout << "--var " << varValue << " ";
- if(stat_opt[0])
- std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << sqrt(varValue) << " ";
+ statfactory::StatFactory stat;
+ vector<double> readBuffer;
+ double varValue;
+ imgReader.readDataBlock(readBuffer, 0, imgReader.nrOfCol()-1, 0, imgReader.nrOfRow()-1, band_opt[0]);
+ stat.setNoDataValues(nodata_opt);
+ stat.meanVar(readBuffer,meanValue,varValue);
+ medianValue=stat.median(readBuffer);
+ stat.minmax(readBuffer,readBuffer.begin(),readBuffer.end(),minValue,maxValue);
+ if(mean_opt[0])
+ std::cout << "--mean " << meanValue << " ";
+ if(median_opt[0])
+ std::cout << "--median " << medianValue << " ";
+ if(stdev_opt[0])
+ std::cout << "--stdDev " << sqrt(varValue) << " ";
+ if(var_opt[0])
+ std::cout << "--var " << varValue << " ";
+ if(stat_opt[0])
+ std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << sqrt(varValue) << " ";
}
if(fstat_opt[0]){//the fast way
- assert(band_opt[iband]<imgReader.nrOfBand());
- GDALProgressFunc pfnProgress;
- void* pProgressData;
- GDALRasterBand* rasterBand;
- rasterBand=imgReader.getRasterBand(band_opt[iband]);
- rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
-
- std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " ";
+ assert(band_opt[iband]<imgReader.nrOfBand());
+ GDALProgressFunc pfnProgress;
+ void* pProgressData;
+ GDALRasterBand* rasterBand;
+ rasterBand=imgReader.getRasterBand(band_opt[iband]);
+ rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
+
+ std::cout << "-min " << minValue << " -max " << maxValue << " --mean " << meanValue << " --stdDev " << stdDev << " ";
}
if(minmax_opt[0]||min_opt[0]||max_opt[0]){
- assert(band_opt[iband]<imgReader.nrOfBand());
-
- if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
- double uli,ulj,lri,lrj;
- imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
- imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
- imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue);
- }
- else{
- imgReader.getMinMax(minValue,maxValue,band_opt[iband]);
- }
- if(minmax_opt[0])
- std::cout << "-min " << minValue << " -max " << maxValue << " ";
- else{
- if(min_opt[0])
- std::cout << "-min " << minValue << " ";
- if(max_opt[0])
- std::cout << "-max " << maxValue << " ";
- }
+ assert(band_opt[iband]<imgReader.nrOfBand());
+
+ if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
+ double uli,ulj,lri,lrj;
+ imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
+ imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
+ imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue);
+ }
+ else{
+ imgReader.getMinMax(minValue,maxValue,band_opt[iband]);
+ }
+ if(minmax_opt[0])
+ std::cout << "-min " << minValue << " -max " << maxValue << " ";
+ else{
+ if(min_opt[0])
+ std::cout << "-min " << minValue << " ";
+ if(max_opt[0])
+ std::cout << "-max " << maxValue << " ";
+ }
}
}
if(histogram_opt[0]){//aggregate results from multiple inputs, but only calculate for first selected band
assert(band_opt[0]<imgReader.nrOfBand());
nbin=(nbin_opt.size())? nbin_opt[0]:0;
-
+
imgReader.getMinMax(minValue,maxValue,band_opt[0]);
if(src_min_opt.size())
minValue=src_min_opt[0];
if(src_max_opt.size())
maxValue=src_max_opt[0];
if(minValue>=maxValue)
- imgReader.getMinMax(minValue,maxValue,band_opt[0]);
+ imgReader.getMinMax(minValue,maxValue,band_opt[0]);
if(verbose_opt[0])
- cout << "number of valid pixels in image: " << imgReader.getNvalid(band_opt[0]) << endl;
+ cout << "number of valid pixels in image: " << imgReader.getNvalid(band_opt[0]) << endl;
nsample+=imgReader.getHistogram(histogramOutput,minValue,maxValue,nbin,band_opt[0],kde_opt[0]);
//only output for last input file
if(ifile==input_opt.size()-1){
- std::cout.precision(10);
- for(int bin=0;bin<nbin;++bin){
- double binValue=0;
- if(nbin==maxValue-minValue+1)
- binValue=minValue+bin;
- else
- binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin;
- std::cout << binValue << " ";
- if(relative_opt[0]||kde_opt[0])
- std::cout << 100.0*static_cast<double>(histogramOutput[bin])/static_cast<double>(nsample) << std::endl;
- else
- std::cout << static_cast<double>(histogramOutput[bin]) << std::endl;
- }
+ std::cout.precision(10);
+ for(int bin=0;bin<nbin;++bin){
+ double binValue=0;
+ if(nbin==maxValue-minValue+1)
+ binValue=minValue+bin;
+ else
+ binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin;
+ std::cout << binValue << " ";
+ if(relative_opt[0]||kde_opt[0])
+ std::cout << 100.0*static_cast<double>(histogramOutput[bin])/static_cast<double>(nsample) << std::endl;
+ else
+ std::cout << static_cast<double>(histogramOutput[bin]) << std::endl;
+ }
}
}
if(histogram2d_opt[0]&&input_opt.size()<2){
@@ -354,80 +354,80 @@ int main(int argc, char *argv[])
imgReader.getMinMax(minX,maxX,band_opt[0]);
imgReader.getMinMax(minY,maxY,band_opt[1]);
if(src_min_opt.size()){
- minX=src_min_opt[0];
- minY=src_min_opt[1];
+ minX=src_min_opt[0];
+ minY=src_min_opt[1];
}
if(src_max_opt.size()){
- maxX=src_max_opt[0];
- maxY=src_max_opt[1];
+ maxX=src_max_opt[0];
+ maxY=src_max_opt[1];
}
nbin=(nbin_opt.size())? nbin_opt[0]:0;
if(nbin<=1){
- std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
- if(minX>=maxX)
- imgReader.getMinMax(minX,maxX,band_opt[0]);
- if(minY>=maxY)
- imgReader.getMinMax(minY,maxY,band_opt[1]);
-
- minValue=(minX<minY)? minX:minY;
- maxValue=(maxX>maxY)? maxX:maxY;
- if(verbose_opt[0])
- std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
- nbin=maxValue-minValue+1;
+ std::cerr << "Warning: number of bins not defined, calculating bins from min and max value" << std::endl;
+ if(minX>=maxX)
+ imgReader.getMinMax(minX,maxX,band_opt[0]);
+ if(minY>=maxY)
+ imgReader.getMinMax(minY,maxY,band_opt[1]);
+
+ minValue=(minX<minY)? minX:minY;
+ maxValue=(maxX>maxY)? maxX:maxY;
+ if(verbose_opt[0])
+ std::cout << "min and max values: " << minValue << ", " << maxValue << std::endl;
+ nbin=maxValue-minValue+1;
}
assert(nbin>1);
double sigma=0;
//kernel density estimation as in http://en.wikipedia.org/wiki/Kernel_density_estimation
if(kde_opt[0]){
- assert(band_opt[0]<imgReader.nrOfBand());
- assert(band_opt[1]<imgReader.nrOfBand());
- GDALProgressFunc pfnProgress;
- void* pProgressData;
- GDALRasterBand* rasterBand;
- double stdDev1=0;
- double stdDev2=0;
- rasterBand=imgReader.getRasterBand(band_opt[0]);
- rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData);
- rasterBand=imgReader.getRasterBand(band_opt[1]);
- rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData);
-
- double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
- if(random_opt[0]>0)
- estimatedSize*=random_opt[0]/100.0;
+ assert(band_opt[0]<imgReader.nrOfBand());
+ assert(band_opt[1]<imgReader.nrOfBand());
+ GDALProgressFunc pfnProgress;
+ void* pProgressData;
+ GDALRasterBand* rasterBand;
+ double stdDev1=0;
+ double stdDev2=0;
+ rasterBand=imgReader.getRasterBand(band_opt[0]);
+ rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData);
+ rasterBand=imgReader.getRasterBand(band_opt[1]);
+ rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData);
+
+ double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
+ if(random_opt[0]>0)
+ estimatedSize*=random_opt[0]/100.0;
sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2);
}
assert(nbin);
if(verbose_opt[0]){
- if(sigma>0)
- std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for bands " << band_opt[0] << " and " << band_opt[1] << std::endl;
- else
- std::cout << "calculating 2d histogram for bands " << band_opt[0] << " and " << band_opt[1] << std::endl;
- std::cout << "nbin: " << nbin << std::endl;
+ if(sigma>0)
+ std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for bands " << band_opt[0] << " and " << band_opt[1] << std::endl;
+ else
+ std::cout << "calculating 2d histogram for bands " << band_opt[0] << " and " << band_opt[1] << std::endl;
+ std::cout << "nbin: " << nbin << std::endl;
}
vector< vector<double> > output;
if(maxX<=minX)
- imgReader.getMinMax(minX,maxX,band_opt[0]);
+ imgReader.getMinMax(minX,maxX,band_opt[0]);
if(maxY<=minY)
- imgReader.getMinMax(minY,maxY,band_opt[1]);
+ imgReader.getMinMax(minY,maxY,band_opt[1]);
if(maxX<=minX){
- std::ostringstream s;
- s<<"Error: could not calculate distribution (minX>=maxX)";
- throw(s.str());
+ std::ostringstream s;
+ s<<"Error: could not calculate distribution (minX>=maxX)";
+ throw(s.str());
}
if(maxY<=minY){
- std::ostringstream s;
- s<<"Error: could not calculate distribution (minY>=maxY)";
- throw(s.str());
+ std::ostringstream s;
+ s<<"Error: could not calculate distribution (minY>=maxY)";
+ throw(s.str());
}
output.resize(nbin);
for(int i=0;i<nbin;++i){
- output[i].resize(nbin);
- for(int j=0;j<nbin;++j)
- output[i][j]=0;
+ output[i].resize(nbin);
+ for(int j=0;j<nbin;++j)
+ output[i][j]=0;
}
int binX=0;
int binY=0;
@@ -437,87 +437,87 @@ int main(int argc, char *argv[])
for(int irow=0;irow<imgReader.nrOfRow();++irow){
if(irow%down_opt[0])
continue;
- imgReader.readData(inputX,irow,band_opt[0]);
- imgReader.readData(inputY,irow,band_opt[1]);
- for(int icol=0;icol<imgReader.nrOfCol();++icol){
+ imgReader.readData(inputX,irow,band_opt[0]);
+ imgReader.readData(inputY,irow,band_opt[1]);
+ for(int icol=0;icol<imgReader.nrOfCol();++icol){
if(icol%down_opt[0])
continue;
- if(random_opt[0]>0){
- double p=static_cast<double>(rand())/(RAND_MAX);
- p*=100.0;
- if(p>random_opt[0])
- continue;//do not select for now, go to next column
- }
- if(imgReader.isNoData(inputX[icol]))
- continue;
- if(imgReader.isNoData(inputY[icol]))
- continue;
- ++nvalid;
- if(inputX[icol]>=maxX)
- binX=nbin-1;
- else if(inputX[icol]<=minX)
- binX=0;
- else
- binX=static_cast<int>(static_cast<double>(inputX[icol]-minX)/(maxX-minX)*nbin);
- if(inputY[icol]>=maxY)
- binY=nbin-1;
- else if(inputY[icol]<=minX)
- binY=0;
- else
- binY=static_cast<int>(static_cast<double>(inputY[icol]-minY)/(maxY-minY)*nbin);
- assert(binX>=0);
- assert(binX<output.size());
- assert(binY>=0);
- assert(binY<output[binX].size());
- if(sigma>0){
- //create kde for Gaussian basis function
- //todo: speed up by calculating first and last bin with non-zero contriubtion...
- for(int ibinX=0;ibinX<nbin;++ibinX){
- double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
- double pdfX=gsl_ran_gaussian_pdf(inputX[icol]-centerX, sigma);
- for(int ibinY=0;ibinY<nbin;++ibinY){
- //calculate \integral_ibinX^(ibinX+1)
- double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
- double pdfY=gsl_ran_gaussian_pdf(inputY[icol]-centerY, sigma);
- output[ibinX][binY]+=pdfX*pdfY;
- }
- }
- }
- else
- ++output[binX][binY];
- }
+ if(random_opt[0]>0){
+ double p=static_cast<double>(rand())/(RAND_MAX);
+ p*=100.0;
+ if(p>random_opt[0])
+ continue;//do not select for now, go to next column
+ }
+ if(imgReader.isNoData(inputX[icol]))
+ continue;
+ if(imgReader.isNoData(inputY[icol]))
+ continue;
+ ++nvalid;
+ if(inputX[icol]>=maxX)
+ binX=nbin-1;
+ else if(inputX[icol]<=minX)
+ binX=0;
+ else
+ binX=static_cast<int>(static_cast<double>(inputX[icol]-minX)/(maxX-minX)*nbin);
+ if(inputY[icol]>=maxY)
+ binY=nbin-1;
+ else if(inputY[icol]<=minX)
+ binY=0;
+ else
+ binY=static_cast<int>(static_cast<double>(inputY[icol]-minY)/(maxY-minY)*nbin);
+ assert(binX>=0);
+ assert(binX<output.size());
+ assert(binY>=0);
+ assert(binY<output[binX].size());
+ if(sigma>0){
+ //create kde for Gaussian basis function
+ //todo: speed up by calculating first and last bin with non-zero contriubtion...
+ for(int ibinX=0;ibinX<nbin;++ibinX){
+ double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
+ double pdfX=gsl_ran_gaussian_pdf(inputX[icol]-centerX, sigma);
+ for(int ibinY=0;ibinY<nbin;++ibinY){
+ //calculate \integral_ibinX^(ibinX+1)
+ double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
+ double pdfY=gsl_ran_gaussian_pdf(inputY[icol]-centerY, sigma);
+ output[ibinX][binY]+=pdfX*pdfY;
+ }
+ }
+ }
+ else
+ ++output[binX][binY];
+ }
}
if(verbose_opt[0])
- cout << "number of valid pixels: " << nvalid << endl;
+ cout << "number of valid pixels: " << nvalid << endl;
for(int binX=0;binX<nbin;++binX){
- cout << endl;
- for(int binY=0;binY<nbin;++binY){
- double binValueX=0;
- if(nbin==maxX-minX+1)
- binValueX=minX+binX;
- else
- binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
- double binValueY=0;
- if(nbin==maxY-minY+1)
- binValueY=minY+binY;
- else
- binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
-
- double value=static_cast<double>(output[binX][binY]);
-
- if(relative_opt[0])
- value*=100.0/nvalid;
-
- cout << binValueX << " " << binValueY << " " << value << std::endl;
- // double value=static_cast<double>(output[binX][binY])/nvalid;
- // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
- }
+ cout << endl;
+ for(int binY=0;binY<nbin;++binY){
+ double binValueX=0;
+ if(nbin==maxX-minX+1)
+ binValueX=minX+binX;
+ else
+ binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
+ double binValueY=0;
+ if(nbin==maxY-minY+1)
+ binValueY=minY+binY;
+ else
+ binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
+
+ double value=static_cast<double>(output[binX][binY]);
+
+ if(relative_opt[0])
+ value*=100.0/nvalid;
+
+ cout << binValueX << " " << binValueY << " " << value << std::endl;
+ // double value=static_cast<double>(output[binX][binY])/nvalid;
+ // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
+ }
}
}
if(reg_opt[0]&&input_opt.size()<2){
if(band_opt.size()<2)
- continue;
+ continue;
imgreg.setDown(down_opt[0]);
imgreg.setThreshold(random_opt[0]);
double c0=0;//offset
@@ -527,7 +527,7 @@ int main(int argc, char *argv[])
}
if(regerr_opt[0]&&input_opt.size()<2){
if(band_opt.size()<2)
- continue;
+ continue;
imgreg.setDown(down_opt[0]);
imgreg.setThreshold(random_opt[0]);
double c0=0;//offset
@@ -537,32 +537,32 @@ int main(int argc, char *argv[])
}
if(rmse_opt[0]&&input_opt.size()<2){
if(band_opt.size()<2)
- continue;
+ continue;
vector<double> xBuffer(imgReader.nrOfCol());
vector<double> yBuffer(imgReader.nrOfCol());
double mse=0;
double nValid=0;
double nPixel=imgReader.nrOfCol()/down_opt[0]*imgReader.nrOfRow()/down_opt[0];
for(int irow;irow<imgReader.nrOfRow();irow+=down_opt[0]){
- imgReader.readData(xBuffer,irow,band_opt[0]);
- imgReader.readData(yBuffer,irow,band_opt[1]);
- for(int icol;icol<imgReader.nrOfCol();icol+=down_opt[0]){
- double xValue=xBuffer[icol];
- double yValue=yBuffer[icol];
- if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
- continue;
- }
- if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
- continue;
- }
- if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[0]||yValue>src_max_opt[0])
- continue;
- ++nValid;
- double e=xValue-yValue;
- if(relative_opt[0])
- e/=yValue;
- mse+=e*e/nPixel;
- }
+ imgReader.readData(xBuffer,irow,band_opt[0]);
+ imgReader.readData(yBuffer,irow,band_opt[1]);
+ for(int icol;icol<imgReader.nrOfCol();icol+=down_opt[0]){
+ double xValue=xBuffer[icol];
+ double yValue=yBuffer[icol];
+ if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
+ continue;
+ }
+ if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
+ continue;
+ }
+ if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[0]||yValue>src_max_opt[0])
+ continue;
+ ++nValid;
+ double e=xValue-yValue;
+ if(relative_opt[0])
+ e/=yValue;
+ mse+=e*e/nPixel;
+ }
}
double correctNorm=nValid;
correctNorm/=nPixel;
@@ -571,7 +571,7 @@ int main(int argc, char *argv[])
}
if(preg_opt[0]&&input_opt.size()<2){
if(band_opt.size()<2)
- continue;
+ continue;
imgreg.setDown(down_opt[0]);
imgreg.setThreshold(random_opt[0]);
double c0=0;//offset
@@ -586,11 +586,11 @@ int main(int argc, char *argv[])
// band_opt.push_back(band_opt[0]);
// if(src_min_opt.size()){
// while(src_min_opt.size()<input_opt.size())
- // src_min_opt.push_back(src_min_opt[0]);
+ // src_min_opt.push_back(src_min_opt[0]);
// }
// if(src_max_opt.size()){
// while(src_max_opt.size()<input_opt.size())
- // src_max_opt.push_back(src_max_opt[0]);
+ // src_max_opt.push_back(src_max_opt[0]);
// }
// ImgReaderGdal imgReader1(input_opt[0]);
// ImgReaderGdal imgReader2(input_opt[1]);
@@ -626,21 +626,21 @@ int main(int argc, char *argv[])
// imgReader1.readData(xBuffer,irow1,band_opt[0]);
// imgReader2.readData(yBuffer,irow2,band_opt[1]);
// for(int icol;icol<imgReader.nrOfCol();icol+=down_opt[0]){
- // icol1=icol;
- // imgReader1.image2geo(icol1,irow1,geoX,geoY);
- // imgReader2.geo2image(geoX,geoY,icol2,irow2);
- // double xValue=xBuffer[icol1];
- // double yValue=yBuffer[icol2];
- // if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
- // continue;
- // }
- // if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[1]||yValue>src_max_opt[1])
- // continue;
- // ++nValid;
- // double e=xValue-yValue;
- // if(relative_opt[0])
- // e/=yValue;
- // mse+=e*e/nPixel;
+ // icol1=icol;
+ // imgReader1.image2geo(icol1,irow1,geoX,geoY);
+ // imgReader2.geo2image(geoX,geoY,icol2,irow2);
+ // double xValue=xBuffer[icol1];
+ // double yValue=yBuffer[icol2];
+ // if(imgReader.isNoData(xValue)||imgReader.isNoData(yValue)){
+ // continue;
+ // }
+ // if(xValue<src_min_opt[0]||xValue>src_max_opt[0]||yValue<src_min_opt[1]||yValue>src_max_opt[1])
+ // continue;
+ // ++nValid;
+ // double e=xValue-yValue;
+ // if(relative_opt[0])
+ // e/=yValue;
+ // mse+=e*e/nPixel;
// }
// }
// double correctNorm=nValid;
@@ -657,11 +657,11 @@ int main(int argc, char *argv[])
band_opt.push_back(band_opt[0]);
if(src_min_opt.size()){
while(src_min_opt.size()<input_opt.size())
- src_min_opt.push_back(src_min_opt[0]);
+ src_min_opt.push_back(src_min_opt[0]);
}
if(src_max_opt.size()){
while(src_max_opt.size()<input_opt.size())
- src_max_opt.push_back(src_max_opt[0]);
+ src_max_opt.push_back(src_max_opt[0]);
}
ImgReaderGdal imgReader1(input_opt[0]);
ImgReaderGdal imgReader2(input_opt[1]);
@@ -698,11 +698,11 @@ int main(int argc, char *argv[])
band_opt.push_back(band_opt[0]);
if(src_min_opt.size()){
while(src_min_opt.size()<input_opt.size())
- src_min_opt.push_back(src_min_opt[0]);
+ src_min_opt.push_back(src_min_opt[0]);
}
if(src_max_opt.size()){
while(src_max_opt.size()<input_opt.size())
- src_max_opt.push_back(src_max_opt[0]);
+ src_max_opt.push_back(src_max_opt[0]);
}
ImgReaderGdal imgReader1(input_opt[0]);
ImgReaderGdal imgReader2(input_opt[1]);
@@ -739,11 +739,11 @@ int main(int argc, char *argv[])
band_opt.push_back(band_opt[0]);
if(src_min_opt.size()){
while(src_min_opt.size()<input_opt.size())
- src_min_opt.push_back(src_min_opt[0]);
+ src_min_opt.push_back(src_min_opt[0]);
}
if(src_max_opt.size()){
while(src_max_opt.size()<input_opt.size())
- src_max_opt.push_back(src_max_opt[0]);
+ src_max_opt.push_back(src_max_opt[0]);
}
ImgReaderGdal imgReader1(input_opt[0]);
ImgReaderGdal imgReader2(input_opt[1]);
@@ -780,11 +780,11 @@ int main(int argc, char *argv[])
band_opt.push_back(band_opt[0]);
if(src_min_opt.size()){
while(src_min_opt.size()<input_opt.size())
- src_min_opt.push_back(src_min_opt[0]);
+ src_min_opt.push_back(src_min_opt[0]);
}
if(src_max_opt.size()){
while(src_max_opt.size()<input_opt.size())
- src_max_opt.push_back(src_max_opt[0]);
+ src_max_opt.push_back(src_max_opt[0]);
}
ImgReaderGdal imgReader1(input_opt[0]);
ImgReaderGdal imgReader2(input_opt[1]);
@@ -817,11 +817,11 @@ int main(int argc, char *argv[])
band_opt.push_back(band_opt[0]);
if(src_min_opt.size()){
while(src_min_opt.size()<input_opt.size())
- src_min_opt.push_back(src_min_opt[0]);
+ src_min_opt.push_back(src_min_opt[0]);
}
if(src_max_opt.size()){
while(src_max_opt.size()<input_opt.size())
- src_max_opt.push_back(src_max_opt[0]);
+ src_max_opt.push_back(src_max_opt[0]);
}
ImgReaderGdal imgReader1(input_opt[0]);
ImgReaderGdal imgReader2(input_opt[1]);
@@ -853,7 +853,7 @@ int main(int argc, char *argv[])
cout << "minY: " << minY << endl;
cout << "maxY: " << maxY << endl;
}
-
+
if(src_min_opt.size()){
minX=src_min_opt[0];
minY=src_min_opt[1];
@@ -869,10 +869,10 @@ int main(int argc, char *argv[])
// imgReader1.getMinMax(minX,maxX,band_opt[0]);
// imgReader2.getMinMax(minY,maxY,band_opt[0]);
if(minX>=maxX)
- imgReader1.getMinMax(minX,maxX,band_opt[0]);
+ imgReader1.getMinMax(minX,maxX,band_opt[0]);
if(minY>=maxY)
- imgReader2.getMinMax(minY,maxY,band_opt[1]);
-
+ imgReader2.getMinMax(minY,maxY,band_opt[1]);
+
minValue=(minX<minY)? minX:minY;
maxValue=(maxX>maxY)? maxX:maxY;
if(verbose_opt[0])
@@ -892,19 +892,19 @@ int main(int argc, char *argv[])
rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev1,pfnProgress,pProgressData);
rasterBand=imgReader2.getRasterBand(band_opt[0]);
rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev2,pfnProgress,pProgressData);
-
+
//todo: think of smarter way how to estimate size (nodata!)
double estimatedSize=1.0*imgReader.getNvalid(band_opt[0])/down_opt[0]/down_opt[0];
if(random_opt[0]>0)
- estimatedSize*=random_opt[0]/100.0;
+ estimatedSize*=random_opt[0]/100.0;
sigma=1.06*sqrt(stdDev1*stdDev2)*pow(estimatedSize,-0.2);
}
assert(nbin);
if(verbose_opt[0]){
if(sigma>0)
- std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl;
+ std::cout << "calculating 2d kernel density estimate with sigma " << sigma << " for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl;
else
- std::cout << "calculating 2d histogram for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl;
+ std::cout << "calculating 2d histogram for datasets " << input_opt[0] << " and " << input_opt[1] << std::endl;
std::cout << "nbin: " << nbin << std::endl;
}
@@ -935,7 +935,7 @@ int main(int argc, char *argv[])
for(int i=0;i<nbin;++i){
output[i].resize(nbin);
for(int j=0;j<nbin;++j)
- output[i][j]=0;
+ output[i][j]=0;
}
int binX=0;
int binY=0;
@@ -950,7 +950,7 @@ int main(int argc, char *argv[])
double irow2=0;
for(int irow=0;irow<imgReader1.nrOfRow();++irow){
if(irow%down_opt[0])
- continue;
+ continue;
irow1=irow;
imgReader1.image2geo(icol1,irow1,geoX,geoY);
imgReader2.geo2image(geoX,geoY,icol2,irow2);
@@ -958,58 +958,58 @@ int main(int argc, char *argv[])
imgReader1.readData(inputX,irow1,band_opt[0]);
imgReader2.readData(inputY,irow2,band_opt[1]);
for(int icol=0;icol<imgReader.nrOfCol();++icol){
- if(icol%down_opt[0])
- continue;
- icol1=icol;
- if(random_opt[0]>0){
- double p=static_cast<double>(rand())/(RAND_MAX);
- p*=100.0;
- if(p>random_opt[0])
- continue;//do not select for now, go to next column
- }
- if(imgReader1.isNoData(inputX[icol]))
- continue;
- imgReader1.image2geo(icol1,irow1,geoX,geoY);
- imgReader2.geo2image(geoX,geoY,icol2,irow2);
- icol2=static_cast<int>(icol2);
- if(imgReader2.isNoData(inputY[icol2]))
- continue;
- // ++nvalid;
- if(inputX[icol1]>=maxX)
- binX=nbin-1;
- else if(inputX[icol]<=minX)
- binX=0;
- else
- binX=static_cast<int>(static_cast<double>(inputX[icol1]-minX)/(maxX-minX)*nbin);
- if(inputY[icol2]>=maxY)
- binY=nbin-1;
- else if(inputY[icol2]<=minY)
- binY=0;
- else
- binY=static_cast<int>(static_cast<double>(inputY[icol2]-minY)/(maxY-minY)*nbin);
- assert(binX>=0);
- assert(binX<output.size());
- assert(binY>=0);
- assert(binY<output[binX].size());
- if(sigma>0){
- //create kde for Gaussian basis function
- //todo: speed up by calculating first and last bin with non-zero contriubtion...
- for(int ibinX=0;ibinX<nbin;++ibinX){
- double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
- double pdfX=gsl_ran_gaussian_pdf(inputX[icol1]-centerX, sigma);
- for(int ibinY=0;ibinY<nbin;++ibinY){
- //calculate \integral_ibinX^(ibinX+1)
- double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
- double pdfY=gsl_ran_gaussian_pdf(inputY[icol2]-centerY, sigma);
- output[ibinX][binY]+=pdfX*pdfY;
- nvalid+=pdfX*pdfY;
- }
- }
- }
- else{
- ++output[binX][binY];
- ++nvalid;
- }
+ if(icol%down_opt[0])
+ continue;
+ icol1=icol;
+ if(random_opt[0]>0){
+ double p=static_cast<double>(rand())/(RAND_MAX);
+ p*=100.0;
+ if(p>random_opt[0])
+ continue;//do not select for now, go to next column
+ }
+ if(imgReader1.isNoData(inputX[icol]))
+ continue;
+ imgReader1.image2geo(icol1,irow1,geoX,geoY);
+ imgReader2.geo2image(geoX,geoY,icol2,irow2);
+ icol2=static_cast<int>(icol2);
+ if(imgReader2.isNoData(inputY[icol2]))
+ continue;
+ // ++nvalid;
+ if(inputX[icol1]>=maxX)
+ binX=nbin-1;
+ else if(inputX[icol]<=minX)
+ binX=0;
+ else
+ binX=static_cast<int>(static_cast<double>(inputX[icol1]-minX)/(maxX-minX)*nbin);
+ if(inputY[icol2]>=maxY)
+ binY=nbin-1;
+ else if(inputY[icol2]<=minY)
+ binY=0;
+ else
+ binY=static_cast<int>(static_cast<double>(inputY[icol2]-minY)/(maxY-minY)*nbin);
+ assert(binX>=0);
+ assert(binX<output.size());
+ assert(binY>=0);
+ assert(binY<output[binX].size());
+ if(sigma>0){
+ //create kde for Gaussian basis function
+ //todo: speed up by calculating first and last bin with non-zero contriubtion...
+ for(int ibinX=0;ibinX<nbin;++ibinX){
+ double centerX=minX+static_cast<double>(maxX-minX)*ibinX/nbin;
+ double pdfX=gsl_ran_gaussian_pdf(inputX[icol1]-centerX, sigma);
+ for(int ibinY=0;ibinY<nbin;++ibinY){
+ //calculate \integral_ibinX^(ibinX+1)
+ double centerY=minY+static_cast<double>(maxY-minY)*ibinY/nbin;
+ double pdfY=gsl_ran_gaussian_pdf(inputY[icol2]-centerY, sigma);
+ output[ibinX][binY]+=pdfX*pdfY;
+ nvalid+=pdfX*pdfY;
+ }
+ }
+ }
+ else{
+ ++output[binX][binY];
+ ++nvalid;
+ }
}
}
if(verbose_opt[0])
@@ -1017,24 +1017,24 @@ int main(int argc, char *argv[])
for(int binX=0;binX<nbin;++binX){
cout << endl;
for(int binY=0;binY<nbin;++binY){
- double binValueX=0;
- if(nbin==maxX-minX+1)
- binValueX=minX+binX;
- else
- binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
- double binValueY=0;
- if(nbin==maxY-minY+1)
- binValueY=minY+binY;
- else
- binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
- double value=static_cast<double>(output[binX][binY]);
-
- if(relative_opt[0]||kde_opt[0])
- value*=100.0/nvalid;
-
- cout << binValueX << " " << binValueY << " " << value << std::endl;
- // double value=static_cast<double>(output[binX][binY])/nvalid;
- // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
+ double binValueX=0;
+ if(nbin==maxX-minX+1)
+ binValueX=minX+binX;
+ else
+ binValueX=minX+static_cast<double>(maxX-minX)*(binX+0.5)/nbin;
+ double binValueY=0;
+ if(nbin==maxY-minY+1)
+ binValueY=minY+binY;
+ else
+ binValueY=minY+static_cast<double>(maxY-minY)*(binY+0.5)/nbin;
+ double value=static_cast<double>(output[binX][binY]);
+
+ if(relative_opt[0]||kde_opt[0])
+ value*=100.0/nvalid;
+
+ cout << binValueX << " " << binValueY << " " << value << std::endl;
+ // double value=static_cast<double>(output[binX][binY])/nvalid;
+ // cout << (maxX-minX)*bin/(nbin-1)+minX << " " << (maxY-minY)*bin/(nbin-1)+minY << " " << value << std::endl;
}
}
imgReader1.close();
@@ -1044,7 +1044,7 @@ int main(int argc, char *argv[])
if(!histogram_opt[0]||histogram2d_opt[0])
std::cout << std::endl;
}
-
+
// int nband=(band_opt.size()) ? band_opt.size() : imgReader.nrOfBand();
// const char* pszMessage;
@@ -1067,77 +1067,77 @@ int main(int argc, char *argv[])
// double iimg,jimg;//image coordinates in img image
// for(icol=0;icol<classReader.nrOfCol();++icol){
// if(icol%down_opt[0])
- // continue;
-
-
- // if(rand_opt[0]>0){
- // gsl_rng* r=stat.getRandomGenerator(time(NULL));
- // //todo: init random number generator using time...
- // if(verbose_opt[0])
- // std::cout << "generating " << rand_opt[0] << " random numbers: " << std::endl;
- // for(unsigned int i=0;i<rand_opt[0];++i)
- // std::cout << i << " " << stat.getRandomValue(r,randdist_opt[0],randa_opt[0],randb_opt[0]) << std::endl;
- // }
-
- // imgreg.setDown(down_opt[0]);
- // imgreg.setThreshold(threshold_opt[0]);
- // double c0=0;//offset
- // double c1=1;//scale
- // double err=uncertNodata_opt[0];//start with high initial value in case we do not have first ob err=imgreg.getRMSE(imgReaderModel1,imgReader,c0,c1,verbose_opt[0]);
-
- // int nband=band_opt.size();
- // if(band_opt[0]<0)
- // nband=imgReader.nrOfBand();
- // for(int iband=0;iband<nband;++iband){
- // unsigned short band_opt[iband]=(band_opt[0]<0)? iband : band_opt[iband];
-
- // if(minmax_opt[0]||min_opt[0]||max_opt[0]){
- // assert(band_opt[iband]<imgReader.nrOfBand());
- // if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
- // double uli,ulj,lri,lrj;
- // imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
- // imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
- // imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue);
- // }
- // else
- // imgReader.getMinMax(minValue,maxValue,band_opt[iband],true);
- // if(minmax_opt[0])
- // std::cout << "-min " << minValue << " -max " << maxValue << " ";
- // else{
- // if(min_opt[0])
- // std::cout << "-min " << minValue << " ";
- // if(max_opt[0])
- // std::cout << "-max " << maxValue << " ";
- // }
- // }
- // }
- // if(relative_opt[0])
- // hist_opt[0]=true;
- // if(hist_opt[0]){
- // assert(band_opt[0]<imgReader.nrOfBand());
- // unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
- // std::vector<unsigned long int> output;
- // minValue=0;
- // maxValue=0;
- // //todo: optimize such that getMinMax is only called once...
- // imgReader.getMinMax(minValue,maxValue,band_opt[0]);
-
- // if(src_min_opt.size())
- // minValue=src_min_opt[0];
- // if(src_max_opt.size())
- // maxValue=src_max_opt[0];
- // unsigned long int nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0]);
- // std::cout.precision(10);
- // for(int bin=0;bin<nbin;++bin){
- // double binValue=0;
- // if(nbin==maxValue-minValue+1)
- // binValue=minValue+bin;
- // else
- // binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin;
- // std::cout << binValue << " ";
- // if(relative_opt[0])
- // std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
- // else
- // std::cout << static_cast<double>(output[bin]) << std::endl;
- // }
- // }
+// continue;
+
+
+// if(rand_opt[0]>0){
+// gsl_rng* r=stat.getRandomGenerator(time(NULL));
+// //todo: init random number generator using time...
+// if(verbose_opt[0])
+// std::cout << "generating " << rand_opt[0] << " random numbers: " << std::endl;
+// for(unsigned int i=0;i<rand_opt[0];++i)
+// std::cout << i << " " << stat.getRandomValue(r,randdist_opt[0],randa_opt[0],randb_opt[0]) << std::endl;
+// }
+
+// imgreg.setDown(down_opt[0]);
+// imgreg.setThreshold(threshold_opt[0]);
+// double c0=0;//offset
+// double c1=1;//scale
+// double err=uncertNodata_opt[0];//start with high initial value in case we do not have first ob err=imgreg.getRMSE(imgReaderModel1,imgReader,c0,c1,verbose_opt[0]);
+
+// int nband=band_opt.size();
+// if(band_opt[0]<0)
+// nband=imgReader.nrOfBand();
+// for(int iband=0;iband<nband;++iband){
+// unsigned short band_opt[iband]=(band_opt[0]<0)? iband : band_opt[iband];
+
+// if(minmax_opt[0]||min_opt[0]||max_opt[0]){
+// assert(band_opt[iband]<imgReader.nrOfBand());
+// if((ulx_opt.size()||uly_opt.size()||lrx_opt.size()||lry_opt.size())&&(imgReader.covers(ulx_opt[0],uly_opt[0],lrx_opt[0],lry_opt[0]))){
+// double uli,ulj,lri,lrj;
+// imgReader.geo2image(ulx_opt[0],uly_opt[0],uli,ulj);
+// imgReader.geo2image(lrx_opt[0],lry_opt[0],lri,lrj);
+// imgReader.getMinMax(static_cast<int>(uli),static_cast<int>(lri),static_cast<int>(ulj),static_cast<int>(lrj),band_opt[iband],minValue,maxValue);
+// }
+// else
+// imgReader.getMinMax(minValue,maxValue,band_opt[iband],true);
+// if(minmax_opt[0])
+// std::cout << "-min " << minValue << " -max " << maxValue << " ";
+// else{
+// if(min_opt[0])
+// std::cout << "-min " << minValue << " ";
+// if(max_opt[0])
+// std::cout << "-max " << maxValue << " ";
+// }
+// }
+// }
+// if(relative_opt[0])
+// hist_opt[0]=true;
+// if(hist_opt[0]){
+// assert(band_opt[0]<imgReader.nrOfBand());
+// unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
+// std::vector<unsigned long int> output;
+// minValue=0;
+// maxValue=0;
+// //todo: optimize such that getMinMax is only called once...
+// imgReader.getMinMax(minValue,maxValue,band_opt[0]);
+
+// if(src_min_opt.size())
+// minValue=src_min_opt[0];
+// if(src_max_opt.size())
+// maxValue=src_max_opt[0];
+// unsigned long int nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0]);
+// std::cout.precision(10);
+// for(int bin=0;bin<nbin;++bin){
+// double binValue=0;
+// if(nbin==maxValue-minValue+1)
+// binValue=minValue+bin;
+// else
+// binValue=minValue+static_cast<double>(maxValue-minValue)*(bin+0.5)/nbin;
+// std::cout << binValue << " ";
+// if(relative_opt[0])
+// std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
+// else
+// std::cout << static_cast<double>(output[bin]) << std::endl;
+// }
+// }
diff --git a/src/apps/pkstatprofile.cc b/src/apps/pkstatprofile.cc
index 841dc3b..a6cac97 100644
--- a/src/apps/pkstatprofile.cc
+++ b/src/apps/pkstatprofile.cc
@@ -209,7 +209,7 @@ int main(int argc,char **argv) {
double gt[6];
input.getGeoTransform(gt);
output.setGeoTransform(gt);
-
+
if(colorTable_opt.size()){
if(colorTable_opt[0]!="none"){
if(verbose_opt[0])
@@ -220,7 +220,7 @@ int main(int argc,char **argv) {
}
else if(input.getColorTable()!=NULL)
output.setColorTable(input.getColorTable());
-
+
if(nodata_opt.size()){
for(int iband=0;iband<output.nrOfBand();++iband)
output.GDALSetNoDataValue(nodata_opt[0],iband);
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index ffc4559..cecce4d 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -31,7 +31,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
/**
Class to read a raster dataset in a format supported by GDAL. Data are cached in memory for a number of rows (if memory>0) before read from file.
- This class inherits from ImgRasterGdal, a general raster class to store e.g., filename, number of columns, rows and bands of the dataset.
+ This class inherits from ImgRasterGdal, a general raster class to store e.g., filename, number of columns, rows and bands of the dataset.
**/
class ImgReaderGdal : public virtual ImgRasterGdal
{
@@ -162,25 +162,45 @@ template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, int mi
**/
template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, int minCol, int maxCol, double row, int band, RESAMPLE resample)
{
+ double eps=0.00001;
std::vector<T> readBuffer_upper;
std::vector<T> readBuffer_lower;
+ double upperRow, lowerRow;
if(buffer.size()!=maxCol-minCol+1)
buffer.resize(maxCol-minCol+1);
- double upperRow=row-0.5;
- upperRow=static_cast<int>(upperRow);
- double lowerRow=row+0.5;
- lowerRow=static_cast<int>(lowerRow);
+ /* double upperRow=row-0.5; */
+ /* upperRow=static_cast<int>(upperRow); */
+ /* double lowerRow=row+0.5; */
+ /* lowerRow=static_cast<int>(lowerRow); */
switch(resample){
case(BILINEAR):
- if(lowerRow>=nrOfRow())
- lowerRow=nrOfRow()-1;
+ /* if(lowerRow>=nrOfRow()) */
+ /* lowerRow=nrOfRow()-1; */
+ upperRow=static_cast<int>(row-0.5+eps);
if(upperRow<0)
upperRow=0;
+ if(upperRow>=nrOfRow())
+ upperRow=nrOfRow()-1;
+ lowerRow=upperRow+1.0;
+ if(lowerRow>=nrOfRow())
+ lowerRow=nrOfRow()-1;
+
readData(readBuffer_upper,minCol,maxCol,static_cast<int>(upperRow),band);
readData(readBuffer_lower,minCol,maxCol,static_cast<int>(lowerRow),band);
//do interpolation in y
for(int icol=0;icol<maxCol-minCol+1;++icol){
- buffer[icol]=(lowerRow-row+0.5)*readBuffer_upper[icol]+(1-lowerRow+row-0.5)*readBuffer_lower[icol];
+ /* buffer[icol]=(lowerRow-row+0.5)*readBuffer_upper[icol]+(1-lowerRow+row-0.5)*readBuffer_lower[icol]; */
+ if(!isNoData(readBuffer_upper[icol])){
+ if(!isNoData(readBuffer_lower[icol])){
+ buffer[icol]=(lowerRow-row+0.5)*readBuffer_upper[icol]+(1-lowerRow+row-0.5)*readBuffer_lower[icol];
+ }
+ else{
+ buffer[icol]=readBuffer_upper[icol];
+ }
+ }
+ else{
+ buffer[icol]=readBuffer_lower[icol];
+ }
}
break;
default:
@@ -211,7 +231,7 @@ template<typename T> void ImgReaderGdal::readDataBlock(Vector2d<T>& buffer2d, in
startit+=maxCol-minCol+1;
}
}
-
+
/**
* @param[out] buffer One dimensional vector representing all pixel values read starting from upper left to lower right.
* @param[in] minCol First column from where to start reading (counting starts from 0)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list