Skip to content
Commits on Source (2)
Expand all Show whitespace changes
Inline Side-by-side
CMakeFiles/cmake.check_cache 0 → 100644
View file @ bd6578cc
# This file is generated by cmake for dependency checking of the CMakeCache.txt file
CMakeLists.txt
View file @ bd6578cc
......@@ -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.3)
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}")
......@@ -135,7 +135,6 @@ if(WIN32)
add_definitions(-D_CRT_NONSTDC_NO_WARNING)
add_definitions(-D_SCL_SECURE_NO_WARNINGS)
endif()
if(CMAKE_CXX_FLAGS MATCHES "/W[0-4]")
string(REGEX REPLACE "/W[0-4]" "/W4"
CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
......@@ -220,7 +219,6 @@ if (BUILD_WITH_LIBLAS)
include_directories(${Boost_INCLUDE_DIRS})
add_definitions("-DBOOST_ALL_NO_LIB")
endif()
# include_directories(${BOOST_INCLUDE_DIR})
# if (MSVC)
# set(BOOST_LIBRARIES -LIBPATH:${BOOST_LIB_PATH} libboost_filesystem-vc100-mt-1_56.lib libboost_system-vc100-mt-1_56.lib)
......@@ -357,7 +355,6 @@ if (PKTOOLS_WITH_UTILITIES)
endif()
###############################################################################
###############################################################################
......@@ -379,7 +376,6 @@ if (PKTOOLS_WITH_UTILITIES)
# if (BUILD_WITH_NLOPT)
# install (TARGETS ${PKNLOPTUTILITIES} DESTINATION bin PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE)
# endif(BUILD_WITH_NLOPT)
endif(PKTOOLS_WITH_UTILITIES)
###############################################################################
......@@ -398,11 +394,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)
......
src/algorithms/Filter.cc
View file @ bd6578cc
......@@ -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,7 +165,7 @@ 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);
......@@ -180,7 +180,7 @@ void filter::Filter::dwtCutFrom(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;
......@@ -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,7 +381,7 @@ 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);
......@@ -418,7 +418,7 @@ void filter::Filter::stat(ImgReaderGdal& input, ImgWriterGdal& output, const std
}
}
try{
output.writeData(lineOutput,GDT_Float64,y);
output.writeData(lineOutput,y);
}
catch(string errorstring){
cerr << errorstring << "in line " << y << endl;
......@@ -443,11 +443,18 @@ 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){
try{
pixelInput=lineInput.selectCol(x);
}
catch(...){
std::cout << "error is caught" << std::endl;
exit(1);
}
int ithreshold=0;//threshold to use for percentiles
try{
for(int imethod=0;imethod<methods.size();++imethod){
switch(getFilterType(methods[imethod])){
case(filter::nvalid):
......@@ -488,14 +495,12 @@ void filter::Filter::stats(ImgReaderGdal& input, ImgWriterGdal& output, const ve
}
}
}
for(int imethod=0;imethod<methods.size();++imethod){
try{
output.writeData(lineOutput[imethod],GDT_Float64,y,imethod);
}
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){
......@@ -526,7 +531,7 @@ void filter::Filter::filter(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;
......
src/apps/pkcomposite.cc
View file @ bd6578cc
......@@ -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"
......@@ -578,9 +579,10 @@ int main(int argc, char *argv[])
// int ncol=static_cast<int>(dcol);
// int nrow=static_cast<int>(drow);
int ncol=ceil((maxLRX-minULX)/dx);
int nrow=ceil((maxULY-minLRY)/dy);
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;
......@@ -701,7 +703,7 @@ int main(int argc, char *argv[])
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);
......@@ -1109,6 +1110,11 @@ int main(int argc, char *argv[])
}
}
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]]){
......@@ -1315,4 +1321,3 @@ int main(int argc, char *argv[])
maskReader.close();
imgWriter.close();
}
src/apps/pkcrop.cc
View file @ bd6578cc
......@@ -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;
......
src/apps/pkextractimg.cc
View file @ bd6578cc
......@@ -90,7 +90,7 @@ using namespace std;
int main(int argc, char *argv[])
{
Optionpk<string> image_opt("i", "input", "Raster input dataset containing band information");
Optionpk<string> sample_opt("s", "sample", "OGR vector dataset with features to be extracted from input data. Output will contain features with input band information included. Sample image can also be GDAL raster dataset.");
Optionpk<string> sample_opt("s", "sample", "Raster dataset with features to be extracted from input data. Output will contain features with input band information included.");
Optionpk<string> output_opt("o", "output", "Output sample dataset");
Optionpk<int> class_opt("c", "class", "Class(es) to extract from input sample image. Leave empty to extract all valid data pixels from sample dataset");
Optionpk<float> threshold_opt("t", "threshold", "Probability threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use a single threshold per vector sample layer. If using raster land cover maps as a sample dataset, you can provide a threshold value for each class (e.g. -t 80 -t 60). Use value 100 to select all pixels for selected class(es)", 100);
......
src/apps/pkextractogr.cc
View file @ bd6578cc
This diff is collapsed.
src/apps/pkkalman.cc
View file @ bd6578cc
......@@ -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);
......@@ -445,6 +451,7 @@ int main(int argc,char **argv) {
cerr << "Error opening file " << model_opt[0] << endl;
}
double modEps=0.00001;
double modRow=0;
double modCol=0;
double lowerCol=0;
......@@ -468,7 +475,7 @@ int main(int argc,char **argv) {
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());
// assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
}
int readModelBand=(model_opt.size()==nmodel)? 0:0;
......@@ -489,16 +496,33 @@ int main(int argc,char **argv) {
continue;
}
}
lowerCol=modCol-0.5;
lowerCol=static_cast<int>(lowerCol);
upperCol=modCol+0.5;
upperCol=static_cast<int>(upperCol);
// 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=(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];
......@@ -577,7 +601,7 @@ int main(int argc,char **argv) {
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());
// assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
if(modelmask_opt.size())
imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
......@@ -593,20 +617,37 @@ int main(int argc,char **argv) {
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());
// 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=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=(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;
......@@ -819,7 +860,7 @@ int main(int argc,char **argv) {
imgUpdaterUncert.image2geo(0,irow,geox,geoy);
if(model_opt.size()==nmodel){
imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
// assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
}
......@@ -827,7 +868,7 @@ int main(int argc,char **argv) {
imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
}
assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
// assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
if(modelmask_opt.size()){
......@@ -883,15 +924,31 @@ int main(int argc,char **argv) {
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=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];
}
// 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);
......@@ -900,15 +957,31 @@ int main(int argc,char **argv) {
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=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=(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())
......@@ -1109,6 +1182,7 @@ int main(int argc,char **argv) {
cerr << "Error opening file " << model_opt[0] << endl;
}
double modEps=0.00001;
double modRow=0;
double modCol=0;
double lowerCol=0;
......@@ -1132,7 +1206,7 @@ int main(int argc,char **argv) {
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());
// assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
}
// imgReaderModel1.readData(estReadBuffer,modRow,0,theResample);
int readModelBand=(model_opt.size()==nmodel)? 0:nmodel-1;
......@@ -1153,15 +1227,31 @@ int main(int argc,char **argv) {
continue;
}
}
lowerCol=modCol-0.5;
lowerCol=static_cast<int>(lowerCol);
upperCol=modCol+0.5;
upperCol=static_cast<int>(upperCol);
// 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=(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];
......@@ -1241,7 +1331,7 @@ int main(int argc,char **argv) {
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());
// assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
imgReaderModel1.readData(estReadBuffer,modRow,readModelBand,theResample);
if(modelmask_opt.size())
imgReaderModel1Mask.readData(modelMaskLineBuffer,modRow,readModelMaskBand);
......@@ -1261,15 +1351,31 @@ int main(int argc,char **argv) {
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=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=(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);
......@@ -1485,7 +1591,7 @@ int main(int argc,char **argv) {
imgUpdaterUncert.image2geo(0,irow,geox,geoy);
if(model_opt.size()==nmodel){
imgReaderModel2.geo2image(geox,geoy,modCol,modRow);
assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
// assert(modRow>=0&&modRow<imgReaderModel2.nrOfRow());
imgReaderModel2.readData(model2LineBuffer,modRow,readModel2Band,theResample);
imgReaderModel1.geo2image(geox,geoy,modCol,modRow);
}
......@@ -1494,7 +1600,7 @@ int main(int argc,char **argv) {
imgReaderModel1.readData(model2LineBuffer,modRow,readModel2Band,theResample);
}
assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
// assert(modRow>=0&&modRow<imgReaderModel1.nrOfRow());
imgReaderModel1.readData(model1LineBuffer,modRow,readModel1Band,theResample);
if(modelmask_opt.size()){
imgReaderModel1Mask.readData(model1MaskLineBuffer,modRow,readModel1MaskBand);
......@@ -1550,15 +1656,31 @@ int main(int argc,char **argv) {
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=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=(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);
......@@ -1568,15 +1690,31 @@ int main(int argc,char **argv) {
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=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];
}
// 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())
......@@ -1869,4 +2007,3 @@ int main(int argc,char **argv) {
if(model_opt.size()<nmodel)
imgReaderModel1.close();
}
src/imageclasses/ImgReaderGdal.h
View file @ bd6578cc
......@@ -162,26 +162,46 @@ 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]; */
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:
readData(buffer,minCol,maxCol,static_cast<int>(row),band);
......
test/data/clipped_shp.prj 0 → 100644
View file @ bd6578cc
PROJCS["IRENET95_Irish_Transverse_Mercator",GEOGCS["GCS_IRENET95",DATUM["D_IRENET95",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",53.5],PARAMETER["central_meridian",-8],PARAMETER["scale_factor",0.99982],PARAMETER["false_easting",600000],PARAMETER["false_northing",750000],UNIT["Meter",1]]
\ No newline at end of file
test/data/clipped_shp.qpj 0 → 100644
View file @ bd6578cc
PROJCS["IRENET95 / Irish Transverse Mercator",GEOGCS["IRENET95",DATUM["IRENET95",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6173"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4173"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",53.5],PARAMETER["central_meridian",-8],PARAMETER["scale_factor",0.99982],PARAMETER["false_easting",600000],PARAMETER["false_northing",750000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2157"]]