[pktools] 173/375: support multiple layers in pkdiff and pkextract
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:11 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit a495db367e35d684e5f5cb40f996c02caeecd382
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Fri Jan 3 17:14:24 2014 +0100
support multiple layers in pkdiff and pkextract
---
ChangeLog | 8 +
src/apps/pkcrop.cc | 96 +-
src/apps/pkdiff.cc | 11 +-
src/apps/pkdsm2shadow.cc | 9 +
src/apps/pkextract.cc | 2967 +++++++++++++++++++-------------------
src/apps/pkfilter.cc | 6 +-
src/apps/pkmosaic.cc | 8 +-
src/imageclasses/ImgReaderGdal.h | 46 +
src/imageclasses/ImgWriterOgr.cc | 176 +--
src/imageclasses/ImgWriterOgr.h | 24 +-
10 files changed, 1739 insertions(+), 1612 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index a30d37e..6232853 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -213,7 +213,15 @@ version 2.5
support sqlite: read training file with small case band names b0, b1, ...
- ImgReaderGdal and ImgWriterGdal
re-design of geotransform
+ - ImgReaderGdal
+ support of scaling and offset for reading raster bands
- pkextract
+ support other ogr file formats than ESRI Shape (e.g, using option -f SQLite)
+ support multiple layers
+ deprecated cbox and rbox
+ - pkdiff
+ support other ogr file formats than ESRI Shape (e.g, using option -f SQLite)
+ support multiple layers
- pkeditogr
support other ogr file formats than ESRI Shape (e.g, using option -f SQLite)
- Filter2d.h
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index b812ae8..45a9ffe 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -49,15 +49,15 @@ int main(int argc, char *argv[])
Optionpk<double> ny_opt("ny", "ny", "image size in y to crop (in meter)");
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<int> band_opt("b", "band", "band index to crop (-1: crop all bands)", -1);
+ Optionpk<int> band_opt("b", "band", "band index to crop (leave empty to retain all bands)");
Optionpk<double> autoscale_opt("as", "autoscale", "scale output to min and max, e.g., --autoscale 0 --autoscale 255");
- Optionpk<double> scale_opt("s", "scale", "output=scale*input+offset", 1);
- Optionpk<double> offset_opt("off", "offset", "output=scale*input+offset", 0);
+ Optionpk<double> scale_opt("s", "scale", "output=scale*input+offset");
+ Optionpk<double> offset_opt("off", "offset", "output=scale*input+offset");
Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
- Optionpk<short> nodata_opt("nodata", "nodata", "Nodata value to put in image if out of bounds.");
+ Optionpk<double> nodata_opt("nodata", "nodata", "Nodata value to put in image if out of bounds.");
Optionpk<string> resample_opt("r", "resampling-method", "Resampling method (near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<string> description_opt("d", "description", "Set image description");
Optionpk<bool> verbose_opt("v", "verbose", "verbose", false);
@@ -111,7 +111,7 @@ int main(int argc, char *argv[])
exit(0);
}
- short nodataValue=nodata_opt.size()? nodata_opt[0] : 0;
+ double nodataValue=nodata_opt.size()? nodata_opt[0] : 0;
RESAMPLE theResample;
if(resample_opt[0]=="near"){
theResample=NEAR;
@@ -149,7 +149,7 @@ int main(int argc, char *argv[])
if(!iimg||imgReader.getDeltaY()<dy)
dy=imgReader.getDeltaY();
}
- if(band_opt[0]>=0)
+ if(band_opt.size())
ncropband+=band_opt.size();
else
ncropband+=imgReader.nrOfBand();
@@ -218,16 +218,20 @@ int main(int argc, char *argv[])
//determine number of output bands
int writeBand=0;//write band
- while(scale_opt.size()<band_opt.size())
- scale_opt.push_back(scale_opt[0]);
- while(offset_opt.size()<band_opt.size())
- offset_opt.push_back(offset_opt[0]);
+ if(scale_opt.size()){
+ while(scale_opt.size()<band_opt.size())
+ scale_opt.push_back(scale_opt[0]);
+ }
+ if(offset_opt.size()){
+ while(offset_opt.size()<band_opt.size())
+ offset_opt.push_back(offset_opt[0]);
+ }
if(autoscale_opt.size()){
assert(autoscale_opt.size()%2==0);
- while(autoscale_opt.size()<band_opt.size()*2){
- autoscale_opt.push_back(autoscale_opt[0]);
- autoscale_opt.push_back(autoscale_opt[1]);
- }
+ // while(autoscale_opt.size()<band_opt.size()*2){
+ // autoscale_opt.push_back(autoscale_opt[0]);
+ // autoscale_opt.push_back(autoscale_opt[1]);
+ // }
}
for(int iimg=0;iimg<input_opt.size();++iimg){
@@ -405,9 +409,9 @@ int main(int argc, char *argv[])
int readncol=endCol-startCol+1;
vector<double> readBuffer(readncol+1);
- int nband=(band_opt[0]<0)?imgReader.nrOfBand():band_opt.size();
+ int nband=(band_opt.size())?band_opt.size() : imgReader.nrOfBand();
for(int iband=0;iband<nband;++iband){
- int readBand=(band_opt[0]<0)?iband:band_opt[iband];
+ int readBand=(band_opt.size()>iband)?band_opt[iband]:iband;
if(verbose_opt[0]){
cout << "extracting band " << readBand << endl;
pfnProgress(progress,pszMessage,pProgressArg);
@@ -418,13 +422,31 @@ int main(int argc, char *argv[])
imgReader.getMinMax(static_cast<int>(startCol),static_cast<int>(endCol),static_cast<int>(startRow),static_cast<int>(endRow),readBand,theMin,theMax);
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);
+ }
+ }
+
double readRow=0;
double readCol=0;
double lowerCol=0;
double upperCol=0;
- for(int irow=0;irow<ncroprow;++irow){
+ for(int irow=0;irow<imgWriter.nrOfRow();++irow){
double x=0;
double y=0;
//convert irow to geo
@@ -443,16 +465,19 @@ int main(int argc, char *argv[])
//readRow=0;
//else if(readRow>=imgReader.nrOfRow())
//readRow=imgReader.nrOfRow()-1;
- for(int ib=0;ib<ncropcol;++ib)
+ for(int ib=0;ib<imgWriter.nrOfCol();++ib)
writeBuffer.push_back(nodataValue);
}
else{
+ if(verbose_opt[0])
+ cout << "reading row: " << readRow << endl;
try{
if(endCol<imgReader.nrOfCol()-1)
imgReader.readData(readBuffer,GDT_Float64,startCol,endCol+1,readRow,readBand,theResample);
else
imgReader.readData(readBuffer,GDT_Float64,startCol,endCol,readRow,readBand,theResample);
- for(int ib=0;ib<ncropcol;++ib){
+ // for(int ib=0;ib<ncropcol;++ib){
+ for(int ib=0;ib<imgWriter.nrOfCol();++ib){
assert(imgWriter.image2geo(ib,irow,x,y));
//lookup corresponding row for irow in this file
imgReader.geo2image(x,y,readCol,readRow);
@@ -493,18 +518,16 @@ int main(int argc, char *argv[])
if(!valid)
writeBuffer.push_back(nodataValue);
else{
- double theScale=1;
- double theOffset=0;
- if(autoscale_opt.size()){
- theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin);
- theOffset=autoscale_opt[0]-theScale*theMin;
- //scale=(ubound-lbound)/(max-min)
- //output=scale*(input-min)+lbound
- }
- else{
- theScale=(scale_opt.size()>1)?scale_opt[iband]:scale_opt[0];
- theOffset=(offset_opt.size()>1)?offset_opt[iband]:offset_opt[0];
- }
+ // double theScale=1;
+ // double theOffset=0;
+ // if(autoscale_opt.size()){
+ // theScale=(autoscale_opt[1]-autoscale_opt[0])/(theMax-theMin);
+ // theOffset=autoscale_opt[0]-theScale*theMin;
+ // }
+ // else{
+ // theScale=(scale_opt.size()>1)?scale_opt[iband]:scale_opt[0];
+ // theOffset=(offset_opt.size()>1)?offset_opt[iband]:offset_opt[0];
+ // }
switch(theResample){
case(BILINEAR):
lowerCol=readCol-0.5;
@@ -515,12 +538,14 @@ int main(int argc, char *argv[])
lowerCol=0;
if(upperCol>=imgReader.nrOfCol())
upperCol=imgReader.nrOfCol()-1;
- writeBuffer.push_back((readCol-0.5-lowerCol)*(readBuffer[upperCol-startCol]*theScale+theOffset)+(1-readCol+0.5+lowerCol)*(readBuffer[lowerCol-startCol]*theScale+theOffset));
+ // writeBuffer.push_back((readCol-0.5-lowerCol)*(readBuffer[upperCol-startCol]*theScale+theOffset)+(1-readCol+0.5+lowerCol)*(readBuffer[lowerCol-startCol]*theScale+theOffset));
+ writeBuffer.push_back((readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol]);
break;
default:
readCol=static_cast<int>(readCol);
readCol-=startCol;//we only start reading from startCol
- writeBuffer.push_back(readBuffer[readCol]*theScale+theOffset);
+ // writeBuffer.push_back(readBuffer[readCol]*theScale+theOffset);
+ writeBuffer.push_back(readBuffer[readCol]);
break;
}
}
@@ -532,7 +557,8 @@ int main(int argc, char *argv[])
exit(2);
}
}
- assert(writeBuffer.size()==ncropcol);
+ if(writeBuffer.size()!=imgWriter.nrOfCol())
+ cout << "writeBuffer.size()=" << writeBuffer.size() << ", imgWriter.nrOfCol()=" << imgWriter.nrOfCol() << endl;
assert(writeBuffer.size()==imgWriter.nrOfCol());
try{
imgWriter.writeData(writeBuffer,GDT_Float64,irow,writeBand);
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index 3a89ec8..d47a40d 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -241,11 +241,11 @@ int main(int argc, char *argv[])
int nlayer=referenceReaderOgr.getDataSource()->GetLayerCount();
for(int ilayer=0;ilayer<nlayer;++ilayer){
progress=0;
- cout << "processing layer " << ilayer << endl;
- if(!verbose_opt[0])
- pfnProgress(progress,pszMessage,pProgressArg);
OGRLayer *readLayer;
readLayer = referenceReaderOgr.getDataSource()->GetLayer(ilayer);
+ cout << "processing layer " << readLayer->GetName() << endl;
+ if(!verbose_opt[0])
+ pfnProgress(progress,pszMessage,pProgressArg);
readLayer->ResetReading();
OGRLayer *writeLayer;
if(output_opt.size()){
@@ -253,11 +253,10 @@ int main(int argc, char *argv[])
cout << "creating output vector file " << output_opt[0] << endl;
// assert(output_opt[0].find(".shp")!=string::npos);
char **papszOptions=NULL;
- string layername=type2string<int>(ilayer);//output_opt[0].substr(0,output_opt[0].find(".shp"));
if(verbose_opt[0])
- cout << "creating layer: " << layername << endl;
+ cout << "creating layer: " << readLayer->GetName() << endl;
// if(ogrWriter.createLayer(layername, referenceReaderOgr.getProjection(ilayer), referenceReaderOgr.getGeometryType(ilayer), papszOptions)==NULL)
- writeLayer=ogrWriter.createLayer(layername, referenceReaderOgr.getProjection(ilayer), wkbPoint, papszOptions);
+ writeLayer=ogrWriter.createLayer(readLayer->GetName(), referenceReaderOgr.getProjection(ilayer), wkbPoint, papszOptions);
assert(writeLayer);
if(verbose_opt[0]){
cout << "created layer" << endl;
diff --git a/src/apps/pkdsm2shadow.cc b/src/apps/pkdsm2shadow.cc
index dc050d6..15ede4e 100644
--- a/src/apps/pkdsm2shadow.cc
+++ b/src/apps/pkdsm2shadow.cc
@@ -41,6 +41,8 @@ int main(int argc,char **argv) {
Optionpk<double> sza_opt("sza", "sza", "Sun zenith angle.");
Optionpk<double> saa_opt("saa", "saa", "Sun azimuth angle (N=0 E=90 S=180 W=270).");
Optionpk<int> flag_opt("f", "flag", "Flag to put in image if pixel shadow", 0);
+ Optionpk<double> scale_opt("s", "scale", "scale used for input dsm: height=scale*input+offset");
+ Optionpk<double> offset_opt("off", "offset", "offset used for input dsm: height=scale*input+offset");
Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
@@ -54,6 +56,8 @@ int main(int argc,char **argv) {
sza_opt.retrieveOption(argc,argv);
saa_opt.retrieveOption(argc,argv);
flag_opt.retrieveOption(argc,argv);
+ scale_opt.retrieveOption(argc,argv);
+ offset_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
@@ -74,6 +78,11 @@ int main(int argc,char **argv) {
assert(input_opt.size());
assert(output_opt.size());
input.open(input_opt[0]);
+ if(scale_opt.size())
+ input.setScale(scale_opt[0]);
+ if(offset_opt.size())
+ input.setOffset(offset_opt[0]);
+
// output.open(output_opt[0],input);
GDALDataType theType=GDT_Unknown;
if(verbose_opt[0])
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index 40bb55b..622a350 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -56,8 +56,8 @@ int main(int argc, char *argv[])
Optionpk<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use multiple threshold values (e.g. -t 80 -t 60) if more classes are to be extracted with random selection. Use value 100 to select all pixels for selected class(es)", 100);
Optionpk<double> min_opt("min", "min", "minimum number of samples to select (0)", 0);
Optionpk<short> boundary_opt("bo", "boundary", "boundary for selecting the sample", 1);
- Optionpk<short> rbox_opt("rb", "rbox", "rectangular boundary box (total width in m) to draw around the selected pixel. Can not combined with class option. Use multiple rbox options for multiple boundary boxes. Use value 0 for no box)", 0);
- Optionpk<short> cbox_opt("cbox", "cbox", "circular boundary (diameter in m) to draw around the selected pixel. Can not combined with class option. Use multiple cbox options for multiple boundary boxes. Use value 0 for no box)", 0);
+ // Optionpk<short> rbox_opt("rb", "rbox", "rectangular boundary box (total width in m) to draw around the selected pixel. Can not combined with class option. Use multiple rbox options for multiple boundary boxes. Use value 0 for no box)", 0);
+ // Optionpk<short> cbox_opt("cbox", "cbox", "circular boundary (diameter in m) to draw around the selected pixel. Can not combined with class option. Use multiple cbox options for multiple boundary boxes. Use value 0 for no box)", 0);
Optionpk<short> disc_opt("circ", "circular", "circular disc kernel boundary", 0);
Optionpk<string> ftype_opt("ft", "ftype", "Field type (only Real or Integer)", "Real");
Optionpk<string> ltype_opt("lt", "ltype", "Label type: In16 or String", "Integer");
@@ -85,8 +85,8 @@ int main(int argc, char *argv[])
threshold_opt.retrieveOption(argc,argv);
min_opt.retrieveOption(argc,argv);
boundary_opt.retrieveOption(argc,argv);
- rbox_opt.retrieveOption(argc,argv);
- cbox_opt.retrieveOption(argc,argv);
+ // rbox_opt.retrieveOption(argc,argv);
+ // cbox_opt.retrieveOption(argc,argv);
disc_opt.retrieveOption(argc,argv);
ftype_opt.retrieveOption(argc,argv);
ltype_opt.retrieveOption(argc,argv);
@@ -242,7 +242,19 @@ int main(int argc, char *argv[])
double progress=0;
srandom(time(NULL));
- if((sample_opt[0].find(".tif"))!=std::string::npos){//raster file
+ bool sampleIsRaster=false;
+ ImgReaderOgr sampleReaderOgr;
+ try{
+ sampleReaderOgr.open(sample_opt[0]);
+ //test
+ // int nlayer=sampleReaderOgr.getDataSource()->GetLayerCount();
+ // cout << "nlayer: " << nlayer << endl;
+ }
+ catch(string errorString){
+ sampleIsRaster=true;
+ }
+
+ if(sampleIsRaster){
if(class_opt.empty()){
std::cout << "Warning: no classes selected, if classes must be extracted, set to -1 for all classes using option -c -1" << std::endl;
ImgReaderGdal classReader;
@@ -746,390 +758,398 @@ int main(int argc, char *argv[])
}
}
else{//vector file
+ if(verbose_opt[0]>1)
+ std::cout << "creating image sample writer " << output_opt[0] << std::endl;
+ ImgWriterOgr ogrWriter;
+ ImgWriterOgr ogrTestWriter;
+ ogrWriter.open(output_opt[0],ogrformat_opt[0]);
+ if(test_opt.size()){
if(verbose_opt[0]>1)
- std::cout << "reading position from shape file " << sample_opt[0] << std::endl;
- ImgReaderOgr sampleReader;
- try{
- sampleReader.open(sample_opt[0]);
- }
- catch(std::string errorString){
- std::cout << errorString << std::endl;
- exit(1);
- }
- if(verbose_opt[0]>1)
- std::cout << "creating image sample writer " << output_opt[0] << std::endl;
- ImgWriterOgr ogrWriter;
- ImgWriterOgr ogrTestWriter;
- ogrWriter.open(output_opt[0],ogrformat_opt[0]);
- if(test_opt.size())
- ogrTestWriter.open(test_opt[0],ogrformat_opt[0]);
- char **papszOptions=NULL;
- ostringstream slayer;
- slayer << "training data";
- std::string layername=slayer.str();
+ std::cout << "creating image test writer " << test_opt[0] << std::endl;
+ ogrTestWriter.open(test_opt[0],ogrformat_opt[0]);
+ }
+
+ //support multiple layers
+ int nlayer=sampleReaderOgr.getDataSource()->GetLayerCount();
+ if(verbose_opt[0])
+ std::cout << "number of layers: " << nlayer << endl;
+
+ for(int ilayer=0;ilayer<nlayer;++ilayer){
+ OGRLayer *readLayer;
+ readLayer = sampleReaderOgr.getDataSource()->GetLayer(ilayer);
+ cout << "processing layer " << readLayer->GetName() << endl;
+ readLayer->ResetReading();
+ OGRLayer *writeLayer;
+ OGRLayer *writeTestLayer;
+
if(polygon_opt[0]){
- if(verbose_opt[0])
- std::cout << "create polygons" << std::endl;
- ogrWriter.createLayer(layername, imgReader.getProjection(), wkbPolygon, papszOptions);
+ if(verbose_opt[0])
+ std::cout << "create polygons" << std::endl;
+ char **papszOptions=NULL;
+ writeLayer=ogrWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPolygon, papszOptions);
if(test_opt.size())
- ogrTestWriter.createLayer(layername, imgReader.getProjection(), wkbPolygon, papszOptions);
+ writeTestLayer=ogrTestWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPolygon, papszOptions);
}
else{
- if(verbose_opt[0])
- std::cout << "create points" << std::endl;
- ogrWriter.createLayer(layername, imgReader.getProjection(), wkbPoint, papszOptions);
- if(test_opt.size())
- ogrTestWriter.createLayer(layername, imgReader.getProjection(), wkbPoint, papszOptions);
+ if(verbose_opt[0])
+ std::cout << "create points in layer " << readLayer->GetName() << std::endl;
+ char **papszOptions=NULL;
+
+ writeLayer=ogrWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPoint, papszOptions);
+ if(test_opt.size()){
+ char **papszOptions=NULL;
+ writeTestLayer=ogrTestWriter.createLayer(readLayer->GetName(), imgReader.getProjection(), wkbPoint, papszOptions);
+ }
+ }
+ if(verbose_opt[0])
+ std::cout << "copy fields" << std::flush << std::endl;
+ ogrWriter.copyFields(sampleReaderOgr,ilayer);
+
+ if(test_opt.size()){
+ if(verbose_opt[0])
+ std::cout << "copy fields test writer" << std::flush << std::endl;
+ ogrTestWriter.copyFields(sampleReaderOgr,ilayer);
}
- ogrWriter.copyFields(sampleReader);
- if(test_opt.size())
- ogrTestWriter.copyFields(sampleReader);
vector<std::string> fieldnames;
- sampleReader.getFields(fieldnames);
- assert(fieldnames.size()==ogrWriter.getFieldCount());
+ if(verbose_opt[0])
+ std::cout << "get fields" << std::flush << std::endl;
+ sampleReaderOgr.getFields(fieldnames);
+ assert(fieldnames.size()==ogrWriter.getFieldCount(ilayer));
map<std::string,double> pointAttributes;
switch(ruleMap[rule_opt[0]]){
- // switch(rule_opt[0]){
+ // switch(rule_opt[0]){
case(rule::proportion):{//proportion for each class
- assert(class_opt.size());
- for(int iclass=0;iclass<class_opt.size();++iclass){
- ostringstream cs;
- cs << class_opt[iclass];
- ogrWriter.createField(cs.str(),fieldType);
- }
- if(keepFeatures_opt[0])
- ogrWriter.createField("origId",OFTInteger);//the fieldId of the original feature
- break;
+ assert(class_opt.size());
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ ostringstream cs;
+ cs << class_opt[iclass];
+ ogrWriter.createField(cs.str(),fieldType,ilayer);
+ }
+ if(keepFeatures_opt[0])
+ ogrWriter.createField("origId",OFTInteger,ilayer);//the fieldId of the original feature
+ break;
}
case(rule::custom):
case(rule::minimum):
case(rule::maximum):
case(rule::maxvote):
- assert(class_opt.size());
- ogrWriter.createField(label_opt[0],fieldType);
- if(test_opt.size())
- ogrTestWriter.createField(label_opt[0],fieldType);
- break;
+ assert(class_opt.size());
+ ogrWriter.createField(label_opt[0],fieldType,ilayer);
+ if(test_opt.size())
+ ogrTestWriter.createField(label_opt[0],fieldType,ilayer);
+ break;
case(rule::point):
case(rule::mean):
case(rule::centroid):
default:{
- if(keepFeatures_opt[0]){
- ogrWriter.createField("origId",OFTInteger);//the fieldId of the original feature
+ if(keepFeatures_opt[0]){
+ ogrWriter.createField("origId",OFTInteger,ilayer);//the fieldId of the original feature
if(test_opt.size())
- ogrTestWriter.createField("origId",OFTInteger);
+ ogrTestWriter.createField("origId",OFTInteger,ilayer);
}
- for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
- for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
- if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
- continue;
- for(int iband=0;iband<nband;++iband){
- int theBand=(band_opt[0]<0)?iband:band_opt[iband];
- ostringstream fs;
- if(theDim>1)
- fs << fieldname_opt[iband] << "_" << windowJ << "_" << windowI;
- else
- fs << fieldname_opt[iband];
- if(verbose_opt[0]>1)
- std::cout << "creating field " << fs.str() << std::endl;
+ for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
+ for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
+ if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
+ continue;
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt[0]<0)?iband:band_opt[iband];
+ ostringstream fs;
+ if(theDim>1)
+ fs << fieldname_opt[iband] << "_" << windowJ << "_" << windowI;
+ else
+ fs << fieldname_opt[iband];
+ if(verbose_opt[0]>1)
+ std::cout << "creating field " << fs.str() << std::endl;
- ogrWriter.createField(fs.str(),fieldType);
+ ogrWriter.createField(fs.str(),fieldType,ilayer);
if(test_opt.size())
- ogrTestWriter.createField(fs.str(),fieldType);
- }
- }
+ ogrTestWriter.createField(fs.str(),fieldType,ilayer);
+ }
+ }
}
- break;
+ break;
}
}
- OGRLayer *readLayer;
- readLayer = sampleReader.getDataSource()->GetLayer(0);
- OGRLayer *writeLayer;
- OGRLayer *writeTestLayer;
- writeLayer=ogrWriter.getDataSource()->GetLayer(0);
- if(test_opt.size())
- writeTestLayer=ogrTestWriter.getDataSource()->GetLayer(0);
- readLayer->ResetReading();
OGRFeature *readFeature;
- int isample=0;
unsigned long int ifeature=0;
- unsigned long int nfeature=sampleReader.getFeatureCount();
- ImgWriterOgr boxWriter;
- if(rbox_opt[0]>0||cbox_opt[0]>0){
- assert(bufferOutput_opt.size());
- assert(test_opt.empty());//not implemented
- if(verbose_opt[0]>1)
- std::cout << "opening box writer " << bufferOutput_opt[0] << std::endl;
- boxWriter.open(bufferOutput_opt[0]);
- std::string layername="buffer";
- boxWriter.createLayer(layername, imgReader.getProjection(), wkbPolygon);
- std::string fieldname="fid";//number of the point
- if(verbose_opt[0]>1)
- std::cout << "creating field " << fieldname << std::endl;
- // ogrWriter.createField(fieldname,OFTInteger);
- boxWriter.createField(fieldname,OFTInteger);
- }
+ unsigned long int nfeature=sampleReaderOgr.getFeatureCount();
+ // ImgWriterOgr boxWriter;
+ // if(rbox_opt[0]>0||cbox_opt[0]>0){
+ // assert(bufferOutput_opt.size());
+ // assert(test_opt.empty());//not implemented
+ // if(verbose_opt[0]>1)
+ // std::cout << "opening box writer " << bufferOutput_opt[0] << std::endl;
+ // boxWriter.open(bufferOutput_opt[0]);
+ // std::string layername="buffer";
+ // boxWriter.createLayer(layername, imgReader.getProjection(), wkbPolygon);
+ // std::string fieldname="fid";//number of the point
+ // if(verbose_opt[0]>1)
+ // std::cout << "creating field " << fieldname << std::endl;
+ // // ogrWriter.createField(fieldname,OFTInteger,ilayer);
+ // boxWriter.createField(fieldname,OFTInteger,ilayer);
+ // }
progress=0;
pfnProgress(progress,pszMessage,pProgressArg);
while( (readFeature = readLayer->GetNextFeature()) != NULL ){
bool writeTest=false;//write this feature to test_opt[0] instead of output_opt
- if(verbose_opt[0]>0)
- std::cout << "reading feature " << readFeature->GetFID() << std::endl;
- if(threshold_opt[0]>0){//percentual value
- double p=static_cast<double>(random())/(RAND_MAX);
- p*=100.0;
- if(p>threshold_opt[0]){
+ if(verbose_opt[0]>0)
+ std::cout << "reading feature " << readFeature->GetFID() << std::endl;
+ if(threshold_opt[0]>0){//percentual value
+ double p=static_cast<double>(random())/(RAND_MAX);
+ p*=100.0;
+ if(p>threshold_opt[0]){
if(test_opt.size())
writeTest=true;
else
continue;//do not select for now, go to next feature
}
- }
- else{//absolute value
- if(ntotalvalid>-threshold_opt[0]){
+ }
+ else{//absolute value
+ if(ntotalvalid>-threshold_opt[0]){
if(test_opt.size())
writeTest=true;
else
continue;//do not select any more pixels, go to next column feature
- }
- }
- if(verbose_opt[0]>0)
- std::cout << "processing feature " << readFeature->GetFID() << std::endl;
- //get x and y from readFeature
- double x,y;
- OGRGeometry *poGeometry;
- poGeometry = readFeature->GetGeometryRef();
- assert(poGeometry!=NULL);
- try{
- if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ){
- assert(class_opt.size()<=1);//class_opt not implemented for point yet
- OGRPoint *poPoint = (OGRPoint *) poGeometry;
- x=poPoint->getX();
- y=poPoint->getY();
- bool valid=true;
- for(int imask=0;imask<mask_opt.size();++imask){
- double colMask,rowMask;//image coordinates in mask image
- if(mask_opt.size()>1){//multiple masks
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
- continue;
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
- continue;
- else{
- maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[imask]=rowMask;
- assert(maskBuffer.size()==maskReader[imask].nrOfBand());
- }
- }
- // char ivalue=0;
- int ivalue=0;
- if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[imask]);
- else//use same invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[0]);
- if(maskBuffer[imask][colMask]==ivalue){
- valid=false;
- break;
- }
- }
- else if(maskReader.size()){
- maskReader[0].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol()){
- continue;
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
- }
+ }
+ }
+ if(verbose_opt[0]>0)
+ std::cout << "processing feature " << readFeature->GetFID() << std::endl;
+ //get x and y from readFeature
+ double x,y;
+ OGRGeometry *poGeometry;
+ poGeometry = readFeature->GetGeometryRef();
+ assert(poGeometry!=NULL);
+ try{
+ if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ){
+ assert(class_opt.size()<=1);//class_opt not implemented for point yet
+ OGRPoint *poPoint = (OGRPoint *) poGeometry;
+ x=poPoint->getX();
+ y=poPoint->getY();
+ bool valid=true;
+ for(int imask=0;imask<mask_opt.size();++imask){
+ double colMask,rowMask;//image coordinates in mask image
+ if(mask_opt.size()>1){//multiple masks
+ maskReader[imask].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
+ continue;
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
+ continue;
+ else{
+ maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[imask]=rowMask;
+ assert(maskBuffer.size()==maskReader[imask].nrOfBand());
+ }
+ }
+ // char ivalue=0;
+ int ivalue=0;
+ if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[imask]);
+ else//use same invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[0]);
+ if(maskBuffer[imask][colMask]==ivalue){
+ valid=false;
+ break;
+ }
+ }
+ else if(maskReader.size()){
+ maskReader[0].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol()){
+ continue;
+ // cerr << colMask << " out of mask col range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
+ }
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow()){
- continue;
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- }
- else{
- maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[0]=rowMask;
- }
- }
- for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
- if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
- valid=false;
- break;
- }
- }
- }
- }
- if(!valid)
- continue;
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow()){
+ continue;
+ // cerr << rowMask << " out of mask row range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
+ }
+ else{
+ maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[0]=rowMask;
+ }
+ }
+ for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
+ if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
+ valid=false;
+ break;
+ }
+ }
+ }
+ }
+ if(!valid)
+ continue;
- double value;
- double i_centre,j_centre;
- if(geo_opt[0])
- imgReader.geo2image(x,y,i_centre,j_centre);
- else{
- i_centre=x;
- j_centre=y;
- }
- //nearest neighbour
- j_centre=static_cast<int>(j_centre);
- i_centre=static_cast<int>(i_centre);
- //check if j_centre is out of bounds
- if(static_cast<int>(j_centre)<0||static_cast<int>(j_centre)>=imgReader.nrOfRow())
- continue;
- //check if i_centre is out of bounds
- if(static_cast<int>(i_centre)<0||static_cast<int>(i_centre)>=imgReader.nrOfCol())
- continue;
+ double value;
+ double i_centre,j_centre;
+ if(geo_opt[0])
+ imgReader.geo2image(x,y,i_centre,j_centre);
+ else{
+ i_centre=x;
+ j_centre=y;
+ }
+ //nearest neighbour
+ j_centre=static_cast<int>(j_centre);
+ i_centre=static_cast<int>(i_centre);
+ //check if j_centre is out of bounds
+ if(static_cast<int>(j_centre)<0||static_cast<int>(j_centre)>=imgReader.nrOfRow())
+ continue;
+ //check if i_centre is out of bounds
+ if(static_cast<int>(i_centre)<0||static_cast<int>(i_centre)>=imgReader.nrOfCol())
+ continue;
- if(rbox_opt[0]){
- assert(test_opt.empty());//not implemented
- vector< vector<OGRPoint*> > points;
- points.resize(rbox_opt.size());
- if(verbose_opt[0]>1)
- std::cout << "creating rectangular box for sample " << isample << ": ";
- for(int ibox=0;ibox<rbox_opt.size();++ibox){
- int npoint=4;
- if(verbose_opt[0]>1)
- std::cout << ibox << " ";
- points[ibox].resize(npoint+1);
- vector<OGRPoint> pbPoint(npoint+1);
- pbPoint[0].setX(x-0.5*rbox_opt[ibox]);
- pbPoint[0].setY(y+0.5*rbox_opt[ibox]);
- points[ibox][0]=&(pbPoint[0]);//start point UL
- points[ibox][4]=&(pbPoint[0]);//end point
- pbPoint[1].setX(x+0.5*rbox_opt[ibox]);
- pbPoint[1].setY(y+0.5*rbox_opt[ibox]);
- points[ibox][1]=&(pbPoint[1]);//UR
- pbPoint[2].setX(x+0.5*rbox_opt[ibox]);
- pbPoint[2].setY(y-0.5*rbox_opt[ibox]);
- points[ibox][2]=&(pbPoint[2]);//LR
- pbPoint[3].setX(x-0.5*rbox_opt[ibox]);
- pbPoint[3].setY(y-0.5*rbox_opt[ibox]);
- points[ibox][3]=&(pbPoint[3]);//LL
- std::string fieldname="fid";//number of the point
- boxWriter.addRing(points[ibox],fieldname,isample);
- // boxWriter.addLineString(points[ibox],fieldname,isample);
- }
- if(verbose_opt[0]>1)
- std::cout << std::endl;
- }
- if(cbox_opt[0]>0){
- vector< vector<OGRPoint*> > points;
- points.resize(cbox_opt.size());
- if(verbose_opt[0]>1)
- std::cout << "creating circular box ";
- for(int ibox=0;ibox<cbox_opt.size();++ibox){
- int npoint=50;
- if(verbose_opt[0]>1)
- std::cout << ibox << " ";
- points[ibox].resize(npoint+1);
- vector<OGRPoint> pbPoint(npoint+1);
- double radius=cbox_opt[ibox]/2.0;
- double alpha=0;
- for(int ipoint=0;ipoint<npoint;++ipoint){
- alpha=ipoint*2.0*PI/static_cast<double>(npoint);
- pbPoint[ipoint].setX(x+radius*cos(alpha));
- pbPoint[ipoint].setY(y+radius*sin(alpha));
- points[ibox][ipoint]=&(pbPoint[ipoint]);
- }
- alpha=0;
- pbPoint[npoint].setX(x+radius*cos(alpha));
- pbPoint[npoint].setY(y+radius*sin(alpha));
- points[ibox][npoint]=&(pbPoint[npoint]);
- std::string fieldname="fid";//number of the point
- boxWriter.addRing(points[ibox],fieldname,isample);
- // boxWriter.addLineString(points[ibox],fieldname,isample);
- }
- if(verbose_opt[0]>1)
- std::cout << std::endl;
- }
+ // if(rbox_opt[0]){
+ // assert(test_opt.empty());//not implemented
+ // vector< vector<OGRPoint*> > points;
+ // points.resize(rbox_opt.size());
+ // if(verbose_opt[0]>1)
+ // std::cout << "creating rectangular box for sample " << isample << ": ";
+ // for(int ibox=0;ibox<rbox_opt.size();++ibox){
+ // int npoint=4;
+ // if(verbose_opt[0]>1)
+ // std::cout << ibox << " ";
+ // points[ibox].resize(npoint+1);
+ // vector<OGRPoint> pbPoint(npoint+1);
+ // pbPoint[0].setX(x-0.5*rbox_opt[ibox]);
+ // pbPoint[0].setY(y+0.5*rbox_opt[ibox]);
+ // points[ibox][0]=&(pbPoint[0]);//start point UL
+ // points[ibox][4]=&(pbPoint[0]);//end point
+ // pbPoint[1].setX(x+0.5*rbox_opt[ibox]);
+ // pbPoint[1].setY(y+0.5*rbox_opt[ibox]);
+ // points[ibox][1]=&(pbPoint[1]);//UR
+ // pbPoint[2].setX(x+0.5*rbox_opt[ibox]);
+ // pbPoint[2].setY(y-0.5*rbox_opt[ibox]);
+ // points[ibox][2]=&(pbPoint[2]);//LR
+ // pbPoint[3].setX(x-0.5*rbox_opt[ibox]);
+ // pbPoint[3].setY(y-0.5*rbox_opt[ibox]);
+ // points[ibox][3]=&(pbPoint[3]);//LL
+ // std::string fieldname="fid";//number of the point
+ // boxWriter.addRing(points[ibox],fieldname,isample);
+ // // boxWriter.addLineString(points[ibox],fieldname,isample);
+ // }
+ // if(verbose_opt[0]>1)
+ // std::cout << std::endl;
+ // }
+ // if(cbox_opt[0]>0){
+ // vector< vector<OGRPoint*> > points;
+ // points.resize(cbox_opt.size());
+ // if(verbose_opt[0]>1)
+ // std::cout << "creating circular box ";
+ // for(int ibox=0;ibox<cbox_opt.size();++ibox){
+ // int npoint=50;
+ // if(verbose_opt[0]>1)
+ // std::cout << ibox << " ";
+ // points[ibox].resize(npoint+1);
+ // vector<OGRPoint> pbPoint(npoint+1);
+ // double radius=cbox_opt[ibox]/2.0;
+ // double alpha=0;
+ // for(int ipoint=0;ipoint<npoint;++ipoint){
+ // alpha=ipoint*2.0*PI/static_cast<double>(npoint);
+ // pbPoint[ipoint].setX(x+radius*cos(alpha));
+ // pbPoint[ipoint].setY(y+radius*sin(alpha));
+ // points[ibox][ipoint]=&(pbPoint[ipoint]);
+ // }
+ // alpha=0;
+ // pbPoint[npoint].setX(x+radius*cos(alpha));
+ // pbPoint[npoint].setY(y+radius*sin(alpha));
+ // points[ibox][npoint]=&(pbPoint[npoint]);
+ // std::string fieldname="fid";//number of the point
+ // boxWriter.addRing(points[ibox],fieldname,isample);
+ // // boxWriter.addLineString(points[ibox],fieldname,isample);
+ // }
+ // if(verbose_opt[0]>1)
+ // std::cout << std::endl;
+ // }
- OGRFeature *writeFeature;
-
+ OGRFeature *writeFeature;
+ if(verbose_opt[0]>1)
+ std::cout << "create feature " << sample_opt[0] << std::endl;
if(writeTest)
writeFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writeFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
- if(verbose_opt[0]>1)
- std::cout << "copying fields from points " << sample_opt[0] << std::endl;
- if(writeFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "copying fields from points " << sample_opt[0] << std::endl;
+ if(writeFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writeFeature->GetFieldCount() << " fields" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writeFeature->GetFieldCount() << " fields" << std::endl;
- vector<double> windowBuffer;
- for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
- for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
- if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
- continue;
- int j=j_centre+windowJ;
- //check if j is out of bounds
- if(static_cast<int>(j)<0||static_cast<int>(j)>=imgReader.nrOfRow())
- continue;
- int i=i_centre+windowI;
- //check if i is out of bounds
- if(static_cast<int>(i)<0||static_cast<int>(i)>=imgReader.nrOfCol())
- continue;
- if(verbose_opt[0]>1)
- std::cout << "reading image value at " << i << "," << j;
- for(int iband=0;iband<nband;++iband){
- int theBand=(band_opt[0]<0)?iband:band_opt[iband];
- imgReader.readData(value,GDT_Float64,i,j,theBand);
- if(verbose_opt[0]>1)
- std::cout << ": " << value << std::endl;
- ostringstream fs;
- if(theDim>1)
- fs << fieldname_opt[iband] << "_" << windowJ << "_" << windowI;
- else
- fs << fieldname_opt[iband];
- if(verbose_opt[0]>1)
- std::cout << "set field " << fs.str() << " to " << value << std::endl;
- switch( fieldType ){
- case OFTInteger:
- writeFeature->SetField(fs.str().c_str(),static_cast<int>(value));
- break;
- case OFTString:
- {
- ostringstream os;
- os << value;
- writeFeature->SetField(fs.str().c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- writeFeature->SetField(fs.str().c_str(),value);
- break;
- case OFTRealList:{
- int fieldIndex=writeFeature->GetFieldIndex(fs.str().c_str());
- int nCount;
- const double *theList;
- theList=writeFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vector<double> vectorList(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[nCount]=value;
- writeFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- break;
- }
- default://not supported
- assert(0);
- break;
- }
- }
- }
- }
- if(keepFeatures_opt[0])
- writeFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(verbose_opt[0]>1)
- std::cout << "creating point feature" << std::endl;
+ vector<double> windowBuffer;
+ for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
+ for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
+ if(disc_opt[0]&&(windowI*windowI+windowJ*windowJ>(theDim/2)*(theDim/2)))
+ continue;
+ int j=j_centre+windowJ;
+ //check if j is out of bounds
+ if(static_cast<int>(j)<0||static_cast<int>(j)>=imgReader.nrOfRow())
+ continue;
+ int i=i_centre+windowI;
+ //check if i is out of bounds
+ if(static_cast<int>(i)<0||static_cast<int>(i)>=imgReader.nrOfCol())
+ continue;
+ if(verbose_opt[0]>1)
+ std::cout << "reading image value at " << i << "," << j;
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt[0]<0)?iband:band_opt[iband];
+ imgReader.readData(value,GDT_Float64,i,j,theBand);
+ if(verbose_opt[0]>1)
+ std::cout << ": " << value << std::endl;
+ ostringstream fs;
+ if(theDim>1)
+ fs << fieldname_opt[iband] << "_" << windowJ << "_" << windowI;
+ else
+ fs << fieldname_opt[iband];
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fs.str() << " to " << value << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ writeFeature->SetField(fs.str().c_str(),static_cast<int>(value));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << value;
+ writeFeature->SetField(fs.str().c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ writeFeature->SetField(fs.str().c_str(),value);
+ break;
+ case OFTRealList:{
+ int fieldIndex=writeFeature->GetFieldIndex(fs.str().c_str());
+ int nCount;
+ const double *theList;
+ theList=writeFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vector<double> vectorList(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[nCount]=value;
+ writeFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ break;
+ }
+ default://not supported
+ assert(0);
+ break;
+ }
+ }
+ }
+ }
+ if(keepFeatures_opt[0])
+ writeFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(verbose_opt[0]>1)
+ std::cout << "creating point feature" << std::endl;
if(writeTest){
if(writeTestLayer->CreateFeature( writeFeature ) != OGRERR_NONE ){
std::string errorString="Failed to create feature in shapefile";
@@ -1142,300 +1162,300 @@ int main(int argc, char *argv[])
throw(errorString);
}
}
- OGRFeature::DestroyFeature( writeFeature );
- ++isample;
- ++ntotalvalid;
+ OGRFeature::DestroyFeature( writeFeature );
+ // ++isample;
+ ++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
- }//if wkbPoint
- else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
+ }//if wkbPoint
+ else if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
- OGRPolygon readPolygon = *((OGRPolygon *) poGeometry);
- OGRPolygon writePolygon;
- OGRLinearRing writeRing;
- OGRPoint writeCentroidPoint;
- OGRFeature *writePolygonFeature;
- OGRFeature *writeCentroidFeature;
+ OGRPolygon readPolygon = *((OGRPolygon *) poGeometry);
+ OGRPolygon writePolygon;
+ OGRLinearRing writeRing;
+ OGRPoint writeCentroidPoint;
+ OGRFeature *writePolygonFeature;
+ OGRFeature *writeCentroidFeature;
- readPolygon.closeRings();
+ readPolygon.closeRings();
- if(verbose_opt[0]>1)
- std::cout << "get centroid point from polygon" << std::endl;
- readPolygon.Centroid(&writeCentroidPoint);
+ if(verbose_opt[0]>1)
+ std::cout << "get centroid point from polygon" << std::endl;
+ readPolygon.Centroid(&writeCentroidPoint);
- double ulx,uly,lrx,lry;
- double uli,ulj,lri,lrj;
- if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
- ulx=writeCentroidPoint.getX();
- uly=writeCentroidPoint.getY();
- lrx=ulx;
- lry=uly;
- }
- else{
- //get envelope
- if(verbose_opt[0])
- std::cout << "reading envelope for polygon " << ifeature << std::endl;
- OGREnvelope* psEnvelope=new OGREnvelope();
- readPolygon.getEnvelope(psEnvelope);
- ulx=psEnvelope->MinX;
- uly=psEnvelope->MaxY;
- lrx=psEnvelope->MaxX;
- lry=psEnvelope->MinY;
- delete psEnvelope;
- }
- if(geo_opt[0]){
- imgReader.geo2image(ulx,uly,uli,ulj);
- imgReader.geo2image(lrx,lry,lri,lrj);
- }
- else{
- uli=ulx;
- ulj=uly;
- lri=lrx;
- lrj=lry;
- }
- //nearest neighbour
- ulj=static_cast<int>(ulj);
- uli=static_cast<int>(uli);
- lrj=static_cast<int>(lrj);
- lri=static_cast<int>(lri);
- //iterate through all pixels
- if(verbose_opt[0]>1)
- std::cout << "bounding box for polygon feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
+ double ulx,uly,lrx,lry;
+ double uli,ulj,lri,lrj;
+ if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
+ ulx=writeCentroidPoint.getX();
+ uly=writeCentroidPoint.getY();
+ lrx=ulx;
+ lry=uly;
+ }
+ else{
+ //get envelope
+ if(verbose_opt[0])
+ std::cout << "reading envelope for polygon " << ifeature << std::endl;
+ OGREnvelope* psEnvelope=new OGREnvelope();
+ readPolygon.getEnvelope(psEnvelope);
+ ulx=psEnvelope->MinX;
+ uly=psEnvelope->MaxY;
+ lrx=psEnvelope->MaxX;
+ lry=psEnvelope->MinY;
+ delete psEnvelope;
+ }
+ if(geo_opt[0]){
+ imgReader.geo2image(ulx,uly,uli,ulj);
+ imgReader.geo2image(lrx,lry,lri,lrj);
+ }
+ else{
+ uli=ulx;
+ ulj=uly;
+ lri=lrx;
+ lrj=lry;
+ }
+ //nearest neighbour
+ ulj=static_cast<int>(ulj);
+ uli=static_cast<int>(uli);
+ lrj=static_cast<int>(lrj);
+ lri=static_cast<int>(lri);
+ //iterate through all pixels
+ if(verbose_opt[0]>1)
+ std::cout << "bounding box for polygon feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
- if(uli<0||lri>=imgReader.nrOfCol()||ulj<0||ulj>=imgReader.nrOfRow())
- continue;
+ if(uli<0||lri>=imgReader.nrOfCol()||ulj<0||ulj>=imgReader.nrOfRow())
+ continue;
- int nPointPolygon=0;
- if(polygon_opt[0]){
+ int nPointPolygon=0;
+ if(polygon_opt[0]){
if(writeTest)
writePolygonFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
- }
- else if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
+ }
+ else if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
if(writeTest)
writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- //previously here
- vector<double> polyValues;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point):
- case(rule::mean):
- case(rule::centroid):
- default:
- polyValues.resize(nband);
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- assert(class_opt.size());
- polyValues.resize(class_opt.size());
- break;
- }
- for(int index=0;index<polyValues.size();++index)
- polyValues[index]=0;
- OGRPoint thePoint;
- for(int j=ulj;j<=lrj;++j){
- for(int i=uli;i<=lri;++i){
- //check if point is on surface
- double x=0;
- double y=0;
- imgReader.image2geo(i,j,x,y);
- thePoint.setX(x);
- thePoint.setY(y);
- if(readPolygon.Contains(&thePoint)){
- bool valid=true;
- for(int imask=0;imask<mask_opt.size();++imask){
- double colMask,rowMask;//image coordinates in mask image
- if(mask_opt.size()>1){//multiple masks
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
- continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[imask].nrOfCol());
- // }
+ //previously here
+ vector<double> polyValues;
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::point):
+ case(rule::mean):
+ case(rule::centroid):
+ default:
+ polyValues.resize(nband);
+ break;
+ case(rule::proportion):
+ case(rule::custom):
+ case(rule::minimum):
+ case(rule::maximum):
+ case(rule::maxvote):
+ assert(class_opt.size());
+ polyValues.resize(class_opt.size());
+ break;
+ }
+ for(int index=0;index<polyValues.size();++index)
+ polyValues[index]=0;
+ OGRPoint thePoint;
+ for(int j=ulj;j<=lrj;++j){
+ for(int i=uli;i<=lri;++i){
+ //check if point is on surface
+ double x=0;
+ double y=0;
+ imgReader.image2geo(i,j,x,y);
+ thePoint.setX(x);
+ thePoint.setY(y);
+ if(readPolygon.Contains(&thePoint)){
+ bool valid=true;
+ for(int imask=0;imask<mask_opt.size();++imask){
+ double colMask,rowMask;//image coordinates in mask image
+ if(mask_opt.size()>1){//multiple masks
+ maskReader[imask].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
+ continue;
+ // {
+ // cerr << colMask << " out of mask col range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[imask].nrOfCol());
+ // }
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
- continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
- else{
- maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[imask]=rowMask;
- assert(maskBuffer.size()==maskReader[imask].nrOfBand());
- }
- }
- // char ivalue=0;
- int ivalue=0;
- if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[imask]);
- else//use same invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[0]);
- if(maskBuffer[imask][colMask]==ivalue){
- valid=false;
- break;
- }
- }
- else if(maskReader.size()){
- maskReader[0].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
- continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
- // }
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
+ continue;
+ // {
+ // cerr << rowMask << " out of mask row range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
+ // }
+ else{
+ maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[imask]=rowMask;
+ assert(maskBuffer.size()==maskReader[imask].nrOfBand());
+ }
+ }
+ // char ivalue=0;
+ int ivalue=0;
+ if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[imask]);
+ else//use same invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[0]);
+ if(maskBuffer[imask][colMask]==ivalue){
+ valid=false;
+ break;
+ }
+ }
+ else if(maskReader.size()){
+ maskReader[0].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
+ continue;
+ // {
+ // cerr << colMask << " out of mask col range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
+ // }
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
- else{
- maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[0]=rowMask;
- }
- }
- for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
- if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
- valid=false;
- break;
- }
- }
- }
- }
- if(!valid)
- continue;
- //check if within raster image
- if(i<0||i>=imgReader.nrOfCol())
- continue;
- if(j<0||j>=imgReader.nrOfRow())
- continue;
- writeRing.addPoint(&thePoint);
- if(verbose_opt[0]>1)
- std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
- ++nPointPolygon;
- OGRFeature *writePointFeature;
- if(!polygon_opt[0]){
- //create feature
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only create point at centroid
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ continue;
+ // {
+ // cerr << rowMask << " out of mask row range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
+ // }
+ else{
+ maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[0]=rowMask;
+ }
+ }
+ for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
+ if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
+ valid=false;
+ break;
+ }
+ }
+ }
+ }
+ if(!valid)
+ continue;
+ //check if within raster image
+ if(i<0||i>=imgReader.nrOfCol())
+ continue;
+ if(j<0||j>=imgReader.nrOfRow())
+ continue;
+ writeRing.addPoint(&thePoint);
+ if(verbose_opt[0]>1)
+ std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
+ ++nPointPolygon;
+ OGRFeature *writePointFeature;
+ if(!polygon_opt[0]){
+ //create feature
+ if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only create point at centroid
if(writeTest)
writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
- if(verbose_opt[0]>1)
- std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
- if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- writePointFeature->SetGeometry(&thePoint);
- OGRGeometry *updateGeometry;
- updateGeometry = writePointFeature->GetGeometryRef();
- OGRPoint *poPoint = (OGRPoint *) updateGeometry;
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
- }
- }
- if(verbose_opt[0]>1)
- std::cout << "reading image value within polygon at position " << i << "," << j;
- for(int iband=0;iband<nband;++iband){
- int theBand=(band_opt[0]<0)?iband:band_opt[iband];
- double value=0;
- imgReader.readData(value,GDT_Float64,i,j,theBand);
- if(verbose_opt[0]>1)
- std::cout << ": " << value << std::endl;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
- int iclass=0;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://in centroid if polygon_opt==true or all values as points if polygon_opt!=true
- case(rule::centroid):
- default:
- polyValues[iband]=value;
- break;
- case(rule::mean)://mean as polygon if polygon_opt==true or as point in centroid if polygon_opt!=true
- polyValues[iband]+=value;
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- for(iclass=0;iclass<class_opt.size();++iclass){
- if(value==class_opt[iclass]){
- assert(polyValues.size()>iclass);
- polyValues[iclass]+=1;//value
- break;
- }
- }
- break;
- }
- }
- else{
- // ostringstream fs;
- // if(imgReader.nrOfBand()==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
- switch( fieldType ){
- case OFTInteger:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),static_cast<int>(value));
- break;
- case OFTString:
- {
- ostringstream os;
- os << value;
- writePointFeature->SetField(fieldname_opt[iband].c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
- break;
- case OFTRealList:{
- int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
- int nCount;
- const double *theList;
- theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vector<double> vectorList(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[nCount]=value;
- writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- break;
- }
- default://not supported
- assert(0);
- break;
- }
- }
- }
- if(!polygon_opt[0]){
- // if(keepFeatures_opt[0])
- // writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only at centroid)
- if(keepFeatures_opt[0])
- writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- //write feature
- if(verbose_opt[0]>1)
- std::cout << "creating point feature" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
+ if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ writePointFeature->SetGeometry(&thePoint);
+ OGRGeometry *updateGeometry;
+ updateGeometry = writePointFeature->GetGeometryRef();
+ OGRPoint *poPoint = (OGRPoint *) updateGeometry;
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
+ }
+ }
+ if(verbose_opt[0]>1)
+ std::cout << "reading image value within polygon at position " << i << "," << j;
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt[0]<0)?iband:band_opt[iband];
+ double value=0;
+ imgReader.readData(value,GDT_Float64,i,j,theBand);
+ if(verbose_opt[0]>1)
+ std::cout << ": " << value << std::endl;
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
+ int iclass=0;
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::point)://in centroid if polygon_opt==true or all values as points if polygon_opt!=true
+ case(rule::centroid):
+ default:
+ polyValues[iband]=value;
+ break;
+ case(rule::mean)://mean as polygon if polygon_opt==true or as point in centroid if polygon_opt!=true
+ polyValues[iband]+=value;
+ break;
+ case(rule::proportion):
+ case(rule::custom):
+ case(rule::minimum):
+ case(rule::maximum):
+ case(rule::maxvote):
+ for(iclass=0;iclass<class_opt.size();++iclass){
+ if(value==class_opt[iclass]){
+ assert(polyValues.size()>iclass);
+ polyValues[iclass]+=1;//value
+ break;
+ }
+ }
+ break;
+ }
+ }
+ else{
+ // ostringstream fs;
+ // if(imgReader.nrOfBand()==1)
+ // fs << fieldname_opt[0];
+ // else
+ // fs << fieldname_opt[0] << theBand;
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),static_cast<int>(value));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << value;
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
+ break;
+ case OFTRealList:{
+ int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
+ int nCount;
+ const double *theList;
+ theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vector<double> vectorList(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[nCount]=value;
+ writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ break;
+ }
+ default://not supported
+ assert(0);
+ break;
+ }
+ }
+ }
+ if(!polygon_opt[0]){
+ // if(keepFeatures_opt[0])
+ // writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only at centroid)
+ if(keepFeatures_opt[0])
+ writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ //write feature
+ if(verbose_opt[0]>1)
+ std::cout << "creating point feature" << std::endl;
if(writeTest){
if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
std::string errorString="Failed to create feature in test shapefile";
@@ -1448,607 +1468,607 @@ int main(int argc, char *argv[])
throw(errorString);
}
}
- //destroy feature
- OGRFeature::DestroyFeature( writePointFeature );
+ //destroy feature
+ OGRFeature::DestroyFeature( writePointFeature );
++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid(2): " << ntotalvalid << std::endl;
}
}
- ++isample;
+ // ++isample;
}
- }
+ }
}
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
- //do not create if no points found within polygon
- if(!nPointPolygon)
- continue;
- //add ring to polygon
- if(polygon_opt[0]){
- writePolygon.addRing(&writeRing);
- writePolygon.closeRings();
- //write geometry of writePolygon
- writePolygonFeature->SetGeometry(&writePolygon);
- if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- if(verbose_opt[0]>1)
- std::cout << "copying new fields write polygon " << sample_opt[0] << std::endl;
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
- //write polygon feature
- }
- else{//write mean value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean)
- //create feature
- if(verbose_opt[0]>1)
- std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
- if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- writeCentroidFeature->SetGeometry(&writeCentroidPoint);
- OGRGeometry *updateGeometry;
- updateGeometry = writeCentroidFeature->GetGeometryRef();
- assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
- }
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://value at each point (or at centroid of polygon if line is set)
- default:{
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- // ostringstream fs;
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
+ //do not create if no points found within polygon
+ if(!nPointPolygon)
+ continue;
+ //add ring to polygon
+ if(polygon_opt[0]){
+ writePolygon.addRing(&writeRing);
+ writePolygon.closeRings();
+ //write geometry of writePolygon
+ writePolygonFeature->SetGeometry(&writePolygon);
+ if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "copying new fields write polygon " << sample_opt[0] << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
+ //write polygon feature
+ }
+ else{//write mean value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean)
+ //create feature
+ if(verbose_opt[0]>1)
+ std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
+ if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ writeCentroidFeature->SetGeometry(&writeCentroidPoint);
+ OGRGeometry *updateGeometry;
+ updateGeometry = writeCentroidFeature->GetGeometryRef();
+ assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
+ }
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::point)://value at each point (or at centroid of polygon if line is set)
+ default:{
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ // ostringstream fs;
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ int theBand=(band_opt[0]<0)?index:band_opt[index];
+ // if(nband==1)
+ // fs << fieldname_opt[0];
+ // else
+ // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }
- default://not supported
- std::cout << "field type not supported yet..." << std::endl;
- break;
- }
- }
- break;
- }
- case(rule::mean):
- case(rule::centroid):{//mean value (written to centroid of polygon if line is not set)
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //test
- if(ruleMap[rule_opt[0]]==rule::centroid)
- assert(nPointPolygon<=1);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- // ostringstream fs;
- theValue/=nPointPolygon;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }
- default://not supported
- std::cout << "field type not supported yet..." << std::endl;
- break;
- }
- }
- break;
- }
- case(rule::proportion):{//proportion classes
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- // stat.sum(polyValues);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
- fs << class_opt[index];
- if(polygon_opt[0])
- writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- }
- break;
- }
- case(rule::custom):{//custom
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- assert(polyValues.size()==2);//11:broadleaved, 12:coniferous
- if(polyValues[0]>=75)//broadleaved
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
- else if(polyValues[1]>=75)//coniferous
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
- else if(polyValues[0]>25&&polyValues[1]>25)//mixed
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
- else{
- if(verbose_opt[0]){
- std::cout << "No valid value in polyValues..." << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- std::cout << theValue << " ";
- }
- std::cout << std::endl;
- }
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
- }
- break;
- }
- case(rule::minimum):{
- //minimum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for min class
- int minClass=stat.max(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]<minClass)
- minClass=class_opt[iclass];
- }
- }
- if(verbose_opt[0]>0)
- std::cout << "minClass: " << minClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),minClass);
- break;
- }
- case(rule::maximum):{
- //maximum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for max class
- int maxClass=stat.min(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]>maxClass)
- maxClass=class_opt[iclass];
- }
- }
- if(verbose_opt[0]>0)
- std::cout << "maxClass: " << maxClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
- break;
- }
- case(rule::maxvote):{
- //maximum votes in polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for class with maximum votes
- int maxClass=stat.min(class_opt);
- vector<double>::iterator maxit;
- maxit=stat.max(polyValues,polyValues.begin(),polyValues.end());
- int maxIndex=distance(polyValues.begin(),maxit);
- maxClass=class_opt[maxIndex];
- if(verbose_opt[0]>0)
- std::cout << "maxClass: " << maxClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
- break;
- }
- }
- if(polygon_opt[0]){
- if(keepFeatures_opt[0])
- writePolygonFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(verbose_opt[0]>1)
- std::cout << "creating polygon feature" << std::endl;
- if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
- std::string errorString="Failed to create polygon feature in shapefile";
- throw(errorString);
- }
- OGRFeature::DestroyFeature( writePolygonFeature );
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << theValue;
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTRealList:{
+ int fieldIndex;
+ int nCount;
+ const double *theList;
+ vector<double> vectorList;
+ if(polygon_opt[0]){
+ fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ else{
+ fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ break;
+ }
+ default://not supported
+ std::cout << "field type not supported yet..." << std::endl;
+ break;
+ }
+ }
+ break;
+ }
+ case(rule::mean):
+ case(rule::centroid):{//mean value (written to centroid of polygon if line is not set)
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //test
+ if(ruleMap[rule_opt[0]]==rule::centroid)
+ assert(nPointPolygon<=1);
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ // ostringstream fs;
+ theValue/=nPointPolygon;
+ int theBand=(band_opt[0]<0)?index:band_opt[index];
+ // if(nband==1)
+ // fs << fieldname_opt[0];
+ // else
+ // fs << fieldname_opt[0] << theBand;
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << theValue;
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTRealList:{
+ int fieldIndex;
+ int nCount;
+ const double *theList;
+ vector<double> vectorList;
+ if(polygon_opt[0]){
+ fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ else{
+ fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ break;
+ }
+ default://not supported
+ std::cout << "field type not supported yet..." << std::endl;
+ break;
+ }
+ }
+ break;
+ }
+ case(rule::proportion):{//proportion classes
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyValues);
+ // stat.sum(polyValues);
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ ostringstream fs;
+ fs << class_opt[index];
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
+ }
+ break;
+ }
+ case(rule::custom):{//custom
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyValues);
+ assert(polyValues.size()==2);//11:broadleaved, 12:coniferous
+ if(polyValues[0]>=75)//broadleaved
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
+ else if(polyValues[1]>=75)//coniferous
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
+ else if(polyValues[0]>25&&polyValues[1]>25)//mixed
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
+ else{
+ if(verbose_opt[0]){
+ std::cout << "No valid value in polyValues..." << std::endl;
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ std::cout << theValue << " ";
+ }
+ std::cout << std::endl;
+ }
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
+ }
+ break;
+ }
+ case(rule::minimum):{
+ //minimum of polygon
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for min class
+ int minClass=stat.max(class_opt);
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(polyValues[iclass]>0){
+ if(verbose_opt[0]>1)
+ std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
+ if(class_opt[iclass]<minClass)
+ minClass=class_opt[iclass];
+ }
+ }
+ if(verbose_opt[0]>0)
+ std::cout << "minClass: " << minClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),minClass);
+ break;
+ }
+ case(rule::maximum):{
+ //maximum of polygon
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for max class
+ int maxClass=stat.min(class_opt);
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(polyValues[iclass]>0){
+ if(verbose_opt[0]>1)
+ std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
+ if(class_opt[iclass]>maxClass)
+ maxClass=class_opt[iclass];
+ }
+ }
+ if(verbose_opt[0]>0)
+ std::cout << "maxClass: " << maxClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
+ break;
+ }
+ case(rule::maxvote):{
+ //maximum votes in polygon
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for class with maximum votes
+ int maxClass=stat.min(class_opt);
+ vector<double>::iterator maxit;
+ maxit=stat.max(polyValues,polyValues.begin(),polyValues.end());
+ int maxIndex=distance(polyValues.begin(),maxit);
+ maxClass=class_opt[maxIndex];
+ if(verbose_opt[0]>0)
+ std::cout << "maxClass: " << maxClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
+ break;
+ }
+ }
+ if(polygon_opt[0]){
+ if(keepFeatures_opt[0])
+ writePolygonFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(verbose_opt[0]>1)
+ std::cout << "creating polygon feature" << std::endl;
+ if(writeLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create polygon feature in shapefile";
+ throw(errorString);
+ }
+ OGRFeature::DestroyFeature( writePolygonFeature );
++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid(1): " << ntotalvalid << std::endl;
- }
- else{
- if(keepFeatures_opt[0])
- writeCentroidFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(verbose_opt[0]>1)
- std::cout << "creating point feature in centroid" << std::endl;
- if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
- std::string errorString="Failed to create point feature in shapefile";
- throw(errorString);
- }
- OGRFeature::DestroyFeature( writeCentroidFeature );
+ }
+ else{
+ if(keepFeatures_opt[0])
+ writeCentroidFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(verbose_opt[0]>1)
+ std::cout << "creating point feature in centroid" << std::endl;
+ if(writeLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
+ std::string errorString="Failed to create point feature in shapefile";
+ throw(errorString);
+ }
+ OGRFeature::DestroyFeature( writeCentroidFeature );
++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
- }
- }
- }
- else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){//todo: try to use virtual OGRGeometry instead of OGRMultiPolygon and OGRPolygon
- OGRMultiPolygon readPolygon = *((OGRMultiPolygon *) poGeometry);
- OGRPolygon writePolygon;
- OGRLinearRing writeRing;
- OGRPoint writeCentroidPoint;
- OGRFeature *writePolygonFeature;
- OGRFeature *writeCentroidFeature;
+ }
+ }
+ }
+ else if(wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon){//todo: try to use virtual OGRGeometry instead of OGRMultiPolygon and OGRPolygon
+ OGRMultiPolygon readPolygon = *((OGRMultiPolygon *) poGeometry);
+ OGRPolygon writePolygon;
+ OGRLinearRing writeRing;
+ OGRPoint writeCentroidPoint;
+ OGRFeature *writePolygonFeature;
+ OGRFeature *writeCentroidFeature;
- readPolygon.closeRings();
+ readPolygon.closeRings();
- if(verbose_opt[0]>1)
- std::cout << "get centroid point from polygon" << std::endl;
- readPolygon.Centroid(&writeCentroidPoint);
+ if(verbose_opt[0]>1)
+ std::cout << "get centroid point from polygon" << std::endl;
+ readPolygon.Centroid(&writeCentroidPoint);
- double ulx,uly,lrx,lry;
- double uli,ulj,lri,lrj;
- if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
- ulx=writeCentroidPoint.getX();
- uly=writeCentroidPoint.getY();
- lrx=ulx;
- lry=uly;
- }
- else{
- //get envelope
- if(verbose_opt[0])
- std::cout << "reading envelope for polygon " << ifeature << std::endl;
- OGREnvelope* psEnvelope=new OGREnvelope();
- readPolygon.getEnvelope(psEnvelope);
- ulx=psEnvelope->MinX;
- uly=psEnvelope->MaxY;
- lrx=psEnvelope->MaxX;
- lry=psEnvelope->MinY;
- delete psEnvelope;
- }
- if(geo_opt[0]){
- imgReader.geo2image(ulx,uly,uli,ulj);
- imgReader.geo2image(lrx,lry,lri,lrj);
- }
- else{
- uli=ulx;
- ulj=uly;
- lri=lrx;
- lrj=lry;
- }
- //nearest neighbour
- ulj=static_cast<int>(ulj);
- uli=static_cast<int>(uli);
- lrj=static_cast<int>(lrj);
- lri=static_cast<int>(lri);
- //iterate through all pixels
- if(verbose_opt[0]>1)
- std::cout << "bounding box for feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
+ double ulx,uly,lrx,lry;
+ double uli,ulj,lri,lrj;
+ if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
+ ulx=writeCentroidPoint.getX();
+ uly=writeCentroidPoint.getY();
+ lrx=ulx;
+ lry=uly;
+ }
+ else{
+ //get envelope
+ if(verbose_opt[0])
+ std::cout << "reading envelope for polygon " << ifeature << std::endl;
+ OGREnvelope* psEnvelope=new OGREnvelope();
+ readPolygon.getEnvelope(psEnvelope);
+ ulx=psEnvelope->MinX;
+ uly=psEnvelope->MaxY;
+ lrx=psEnvelope->MaxX;
+ lry=psEnvelope->MinY;
+ delete psEnvelope;
+ }
+ if(geo_opt[0]){
+ imgReader.geo2image(ulx,uly,uli,ulj);
+ imgReader.geo2image(lrx,lry,lri,lrj);
+ }
+ else{
+ uli=ulx;
+ ulj=uly;
+ lri=lrx;
+ lrj=lry;
+ }
+ //nearest neighbour
+ ulj=static_cast<int>(ulj);
+ uli=static_cast<int>(uli);
+ lrj=static_cast<int>(lrj);
+ lri=static_cast<int>(lri);
+ //iterate through all pixels
+ if(verbose_opt[0]>1)
+ std::cout << "bounding box for feature " << ifeature << ": " << uli << " " << ulj << " " << lri << " " << lrj << std::endl;
- if(uli<0||lri>=imgReader.nrOfCol()||ulj<0||ulj>=imgReader.nrOfRow())
- continue;
+ if(uli<0||lri>=imgReader.nrOfCol()||ulj<0||ulj>=imgReader.nrOfRow())
+ continue;
- int nPointPolygon=0;
- if(polygon_opt[0]){
+ int nPointPolygon=0;
+ if(polygon_opt[0]){
if(writeTest)
writePolygonFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writePolygonFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
- }
- else if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
+ }
+ else if(ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
if(writeTest)
writeCentroidFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writeCentroidFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
}
- //previously here
- vector<double> polyValues;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point):
- case(rule::mean):
- case(rule::centroid):
- default:
- polyValues.resize(nband);
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- assert(class_opt.size());
- polyValues.resize(class_opt.size());
- break;
- }
- for(int index=0;index<polyValues.size();++index)
- polyValues[index]=0;
- OGRPoint thePoint;
- for(int j=ulj;j<=lrj;++j){
- for(int i=uli;i<=lri;++i){
- //check if point is on surface
- double x=0;
- double y=0;
- imgReader.image2geo(i,j,x,y);
- thePoint.setX(x);
- thePoint.setY(y);
- if(readPolygon.Contains(&thePoint)){
- bool valid=true;
- for(int imask=0;imask<mask_opt.size();++imask){
- double colMask,rowMask;//image coordinates in mask image
- if(mask_opt.size()>1){//multiple masks
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
- continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[imask].nrOfCol());
- // }
+ //previously here
+ vector<double> polyValues;
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::point):
+ case(rule::mean):
+ case(rule::centroid):
+ default:
+ polyValues.resize(nband);
+ break;
+ case(rule::proportion):
+ case(rule::custom):
+ case(rule::minimum):
+ case(rule::maximum):
+ case(rule::maxvote):
+ assert(class_opt.size());
+ polyValues.resize(class_opt.size());
+ break;
+ }
+ for(int index=0;index<polyValues.size();++index)
+ polyValues[index]=0;
+ OGRPoint thePoint;
+ for(int j=ulj;j<=lrj;++j){
+ for(int i=uli;i<=lri;++i){
+ //check if point is on surface
+ double x=0;
+ double y=0;
+ imgReader.image2geo(i,j,x,y);
+ thePoint.setX(x);
+ thePoint.setY(y);
+ if(readPolygon.Contains(&thePoint)){
+ bool valid=true;
+ for(int imask=0;imask<mask_opt.size();++imask){
+ double colMask,rowMask;//image coordinates in mask image
+ if(mask_opt.size()>1){//multiple masks
+ maskReader[imask].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[imask].nrOfCol())
+ continue;
+ // {
+ // cerr << colMask << " out of mask col range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[imask].nrOfCol());
+ // }
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
- continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
- else{
- maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[imask]=rowMask;
- assert(maskBuffer.size()==maskReader[imask].nrOfBand());
- }
- }
- // char ivalue=0;
- int ivalue=0;
- if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[imask]);
- else//use same invalid value for each mask
- ivalue=static_cast<int>(msknodata_opt[0]);
- if(maskBuffer[imask][colMask]==ivalue){
- valid=false;
- break;
- }
- }
- else if(maskReader.size()){
- maskReader[0].geo2image(x,y,colMask,rowMask);
- //nearest neighbour
- rowMask=static_cast<int>(rowMask);
- colMask=static_cast<int>(colMask);
- if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
- continue;
- // {
- // cerr << colMask << " out of mask col range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
- // }
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[imask])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[imask].nrOfRow())
+ continue;
+ // {
+ // cerr << rowMask << " out of mask row range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
+ // }
+ else{
+ maskReader[imask].readData(maskBuffer[imask],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[imask]=rowMask;
+ assert(maskBuffer.size()==maskReader[imask].nrOfBand());
+ }
+ }
+ // char ivalue=0;
+ int ivalue=0;
+ if(mask_opt.size()==msknodata_opt.size())//one invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[imask]);
+ else//use same invalid value for each mask
+ ivalue=static_cast<int>(msknodata_opt[0]);
+ if(maskBuffer[imask][colMask]==ivalue){
+ valid=false;
+ break;
+ }
+ }
+ else if(maskReader.size()){
+ maskReader[0].geo2image(x,y,colMask,rowMask);
+ //nearest neighbour
+ rowMask=static_cast<int>(rowMask);
+ colMask=static_cast<int>(colMask);
+ if(static_cast<int>(colMask)<0||static_cast<int>(colMask)>=maskReader[0].nrOfCol())
+ continue;
+ // {
+ // cerr << colMask << " out of mask col range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(colMask)>=0&&static_cast<int>(colMask)<maskReader[0].nrOfCol());
+ // }
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- // {
- // cerr << rowMask << " out of mask row range!" << std::endl;
- // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
- // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
- // }
- else{
- maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
- oldmaskrow[0]=rowMask;
- }
- }
- for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
- if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
- valid=false;
- break;
- }
- }
- }
- }
- if(!valid)
- continue;
- //check if within raster image
- if(i<0||i>=imgReader.nrOfCol())
- continue;
- if(j<0||j>=imgReader.nrOfRow())
- continue;
- writeRing.addPoint(&thePoint);
- if(verbose_opt[0]>1)
- std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
- ++nPointPolygon;
- OGRFeature *writePointFeature;
- if(!polygon_opt[0]){
- //create feature
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only create point at centroid)
+ if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ continue;
+ // {
+ // cerr << rowMask << " out of mask row range!" << std::endl;
+ // cerr << x << " " << y << " " << colMask << " " << rowMask << std::endl;
+ // assert(static_cast<int>(rowMask)>=0&&static_cast<int>(rowMask)<imgReader.nrOfRow());
+ // }
+ else{
+ maskReader[0].readData(maskBuffer[0],GDT_Int32,static_cast<int>(rowMask));
+ oldmaskrow[0]=rowMask;
+ }
+ }
+ for(int ivalue=0;ivalue<msknodata_opt.size();++ivalue){
+ if(maskBuffer[0][colMask]==static_cast<int>(msknodata_opt[ivalue])){
+ valid=false;
+ break;
+ }
+ }
+ }
+ }
+ if(!valid)
+ continue;
+ //check if within raster image
+ if(i<0||i>=imgReader.nrOfCol())
+ continue;
+ if(j<0||j>=imgReader.nrOfRow())
+ continue;
+ writeRing.addPoint(&thePoint);
+ if(verbose_opt[0]>1)
+ std::cout << "point is on surface:" << thePoint.getX() << "," << thePoint.getY() << std::endl;
+ ++nPointPolygon;
+ OGRFeature *writePointFeature;
+ if(!polygon_opt[0]){
+ //create feature
+ if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only create point at centroid)
if(writeTest)
writePointFeature = OGRFeature::CreateFeature(writeTestLayer->GetLayerDefn());
else
writePointFeature = OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
- if(verbose_opt[0]>1)
- std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
- if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- writePointFeature->SetGeometry(&thePoint);
- OGRGeometry *updateGeometry;
- updateGeometry = writePointFeature->GetGeometryRef();
- OGRPoint *poPoint = (OGRPoint *) updateGeometry;
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
- }
- }
- if(verbose_opt[0]>1)
- std::cout << "reading image value withinin polygon at position " << i << "," << j;
- for(int iband=0;iband<nband;++iband){
- int theBand=(band_opt[0]<0)?iband:band_opt[iband];
- double value=0;
- imgReader.readData(value,GDT_Float64,i,j,theBand);
- if(verbose_opt[0]>1)
- std::cout << ": " << value << std::endl;
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
- int iclass=0;
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://in centroid if polygon_opt==true or all values as points if polygon_opt!=true
- case(rule::centroid):
- default:
- polyValues[iband]=value;
- break;
- case(rule::mean)://mean as polygon if polygon_opt==true or as point in centroid if polygon_opt!=true
- polyValues[iband]+=value;
- break;
- case(rule::proportion):
- case(rule::custom):
- case(rule::minimum):
- case(rule::maximum):
- case(rule::maxvote):
- for(iclass=0;iclass<class_opt.size();++iclass){
- if(value==class_opt[iclass]){
- assert(polyValues.size()>iclass);
- polyValues[iclass]+=1;//value
- break;
- }
- }
- break;
- }
- }
- else{
- ostringstream fs;
- // if(imgReader.nrOfBand()==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
- switch( fieldType ){
- case OFTInteger:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),static_cast<int>(value));
- break;
- case OFTString:
- {
- ostringstream os;
- os << value;
- writePointFeature->SetField(fieldname_opt[iband].c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
- break;
- case OFTRealList:{
- int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
- int nCount;
- const double *theList;
- theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vector<double> vectorList(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[nCount]=value;
- writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- break;
- }
- default://not supported
- assert(0);
- break;
- }
- }
- }
- if(!polygon_opt[0]){
- if(keepFeatures_opt[0])
- writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only at centroid)
- //write feature
- if(verbose_opt[0]>1)
- std::cout << "creating point feature" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
+ if(writePointFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ writePointFeature->SetGeometry(&thePoint);
+ OGRGeometry *updateGeometry;
+ updateGeometry = writePointFeature->GetGeometryRef();
+ OGRPoint *poPoint = (OGRPoint *) updateGeometry;
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writePointFeature->GetFieldCount() << " fields" << std::endl;
+ }
+ }
+ if(verbose_opt[0]>1)
+ std::cout << "reading image value withinin polygon at position " << i << "," << j;
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt[0]<0)?iband:band_opt[iband];
+ double value=0;
+ imgReader.readData(value,GDT_Float64,i,j,theBand);
+ if(verbose_opt[0]>1)
+ std::cout << ": " << value << std::endl;
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
+ int iclass=0;
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::point)://in centroid if polygon_opt==true or all values as points if polygon_opt!=true
+ case(rule::centroid):
+ default:
+ polyValues[iband]=value;
+ break;
+ case(rule::mean)://mean as polygon if polygon_opt==true or as point in centroid if polygon_opt!=true
+ polyValues[iband]+=value;
+ break;
+ case(rule::proportion):
+ case(rule::custom):
+ case(rule::minimum):
+ case(rule::maximum):
+ case(rule::maxvote):
+ for(iclass=0;iclass<class_opt.size();++iclass){
+ if(value==class_opt[iclass]){
+ assert(polyValues.size()>iclass);
+ polyValues[iclass]+=1;//value
+ break;
+ }
+ }
+ break;
+ }
+ }
+ else{
+ ostringstream fs;
+ // if(imgReader.nrOfBand()==1)
+ // fs << fieldname_opt[0];
+ // else
+ // fs << fieldname_opt[0] << theBand;
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[iband] << " to " << value << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),static_cast<int>(value));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << value;
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ writePointFeature->SetField(fieldname_opt[iband].c_str(),value);
+ break;
+ case OFTRealList:{
+ int fieldIndex=writePointFeature->GetFieldIndex(fieldname_opt[iband].c_str());
+ int nCount;
+ const double *theList;
+ theList=writePointFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vector<double> vectorList(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[nCount]=value;
+ writePointFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ break;
+ }
+ default://not supported
+ assert(0);
+ break;
+ }
+ }
+ }
+ if(!polygon_opt[0]){
+ if(keepFeatures_opt[0])
+ writePointFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(ruleMap[rule_opt[0]]!=rule::mean&&ruleMap[rule_opt[0]]!=rule::centroid){//do not create in case of mean value (only at centroid)
+ //write feature
+ if(verbose_opt[0]>1)
+ std::cout << "creating point feature" << std::endl;
if(writeTest){
if(writeTestLayer->CreateFeature( writePointFeature ) != OGRERR_NONE ){
std::string errorString="Failed to create feature in shapefile";
@@ -2061,280 +2081,280 @@ int main(int argc, char *argv[])
throw(errorString);
}
}
- //destroy feature
- OGRFeature::DestroyFeature( writePointFeature );
- }
- }
- ++isample;
- ++ntotalvalid;
+ //destroy feature
+ OGRFeature::DestroyFeature( writePointFeature );
+ }
+ }
+ // ++isample;
+ ++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
- }
- }
- }
- if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
- //add ring to polygon
- if(polygon_opt[0]){
- writePolygon.addRing(&writeRing);
- writePolygon.closeRings();
- //write geometry of writePolygon
- writePolygonFeature->SetGeometry(&writePolygon);
- if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- if(verbose_opt[0]>1)
- std::cout << "copying new fields write polygon " << sample_opt[0] << std::endl;
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
- //write polygon feature
- }
- else{//write mean value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean)
- //create feature
- if(verbose_opt[0]>1)
- std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
- if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
- cerr << "writing feature failed" << std::endl;
- writeCentroidFeature->SetGeometry(&writeCentroidPoint);
- OGRGeometry *updateGeometry;
- updateGeometry = writeCentroidFeature->GetGeometryRef();
- assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
- if(verbose_opt[0]>1)
- std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
- }
- switch(ruleMap[rule_opt[0]]){
- case(rule::point)://value at each point (or at centroid of polygon if line is set)
- default:{
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }//case OFTRealList
- }//switch(fieldType)
- }//for index
- break;
- }//case 0 and default
- case(rule::mean):
- case(rule::centroid):{//mean value (written to centroid of polygon if line is not set
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //test
- if(ruleMap[rule_opt[0]]==rule::centroid)
- assert(nPointPolygon<=1);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
- theValue/=nPointPolygon;
- int theBand=(band_opt[0]<0)?index:band_opt[index];
- // if(nband==1)
- // fs << fieldname_opt[0];
- // else
- // fs << fieldname_opt[0] << theBand;
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
- switch( fieldType ){
- case OFTInteger:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
- break;
- case OFTString:
- {
- ostringstream os;
- os << theValue;
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
- break;
- }
- case OFTReal:
- if(polygon_opt[0])
- writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
- else
- writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
- break;
- case OFTRealList:{
- int fieldIndex;
- int nCount;
- const double *theList;
- vector<double> vectorList;
- if(polygon_opt[0]){
- fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- else{
- fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
- theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
- vectorList.resize(nCount+1);
- for(int index=0;index<nCount;++index)
- vectorList[index]=theList[index];
- vectorList[nCount]=theValue;
- writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
- }
- break;
- }
- }
- }
- break;
- }
- case(rule::proportion):{//proportion classes
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- // stat.sum(polyValues);
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- ostringstream fs;
- fs << class_opt[index];
- if(polygon_opt[0])
- writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- else
- writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
- }
- break;
- }
- case(rule::custom):{//custom
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- stat.normalize_pct(polyValues);
- assert(polyValues.size()==2);//11:broadleaved, 12:coniferous
- if(polyValues[0]>=75)//broadleaved
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
- else if(polyValues[1]>=75)//coniferous
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
- else if(polyValues[0]>25&&polyValues[1]>25)//mixed
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
- else{
- if(verbose_opt[0]){
- std::cout << "No valid value in polyValues..." << std::endl;
- for(int index=0;index<polyValues.size();++index){
- double theValue=polyValues[index];
- std::cout << theValue << " ";
- }
- std::cout << std::endl;
- }
- writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
- }
- break;
- }
- case(rule::minimum):{//minimum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for min class
- int minClass=stat.max(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]<minClass)
- minClass=class_opt[iclass];
- }
- }
- if(verbose_opt[0]>0)
- std::cout << "minClass: " << minClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),minClass);
- break;
- }
- case(rule::maximum):{//maximum of polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for max class
- int maxClass=stat.min(class_opt);
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(polyValues[iclass]>0){
- if(verbose_opt[0]>1)
- std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
- if(class_opt[iclass]>maxClass)
- maxClass=class_opt[iclass];
- }
- }
- if(verbose_opt[0]>0)
- std::cout << "maxClass: " << maxClass << std::endl;
- writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
- break;
- }
- case(rule::maxvote):{//maximum votes in polygon
- assert(polygon_opt[0]);//not implemented for points
- if(verbose_opt[0])
- std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
- //search for max votes
- int maxClass=stat.min(class_opt);
- vector<double>::iterator maxit;
- maxit=stat.max(polyValues,polyValues.begin(),polyValues.end());
- int maxIndex=distance(polyValues.begin(),maxit);
- maxClass=class_opt[maxIndex];
- }
- }
- if(polygon_opt[0]){
- if(keepFeatures_opt[0])
- writePolygonFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(verbose_opt[0]>1)
- std::cout << "creating polygon feature" << std::endl;
+ }
+ }
+ }
+ if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean||ruleMap[rule_opt[0]]==rule::centroid){
+ //add ring to polygon
+ if(polygon_opt[0]){
+ writePolygon.addRing(&writeRing);
+ writePolygon.closeRings();
+ //write geometry of writePolygon
+ writePolygonFeature->SetGeometry(&writePolygon);
+ if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "copying new fields write polygon " << sample_opt[0] << std::endl;
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writePolygonFeature->GetFieldCount() << " fields" << std::endl;
+ //write polygon feature
+ }
+ else{//write mean value of polygon to centroid point (ruleMap[rule_opt[0]]==rule::mean)
+ //create feature
+ if(verbose_opt[0]>1)
+ std::cout << "copying fields from polygons " << sample_opt[0] << std::endl;
+ if(writeCentroidFeature->SetFrom(readFeature)!= OGRERR_NONE)
+ cerr << "writing feature failed" << std::endl;
+ writeCentroidFeature->SetGeometry(&writeCentroidPoint);
+ OGRGeometry *updateGeometry;
+ updateGeometry = writeCentroidFeature->GetGeometryRef();
+ assert(wkbFlatten(updateGeometry->getGeometryType()) == wkbPoint );
+ if(verbose_opt[0]>1)
+ std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
+ }
+ switch(ruleMap[rule_opt[0]]){
+ case(rule::point)://value at each point (or at centroid of polygon if line is set)
+ default:{
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ ostringstream fs;
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ int theBand=(band_opt[0]<0)?index:band_opt[index];
+ // if(nband==1)
+ // fs << fieldname_opt[0];
+ // else
+ // fs << fieldname_opt[0] << theBand;
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << theValue;
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTRealList:{
+ int fieldIndex;
+ int nCount;
+ const double *theList;
+ vector<double> vectorList;
+ if(polygon_opt[0]){
+ fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ else{
+ fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ break;
+ }//case OFTRealList
+ }//switch(fieldType)
+ }//for index
+ break;
+ }//case 0 and default
+ case(rule::mean):
+ case(rule::centroid):{//mean value (written to centroid of polygon if line is not set
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //test
+ if(ruleMap[rule_opt[0]]==rule::centroid)
+ assert(nPointPolygon<=1);
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ ostringstream fs;
+ theValue/=nPointPolygon;
+ int theBand=(band_opt[0]<0)?index:band_opt[index];
+ // if(nband==1)
+ // fs << fieldname_opt[0];
+ // else
+ // fs << fieldname_opt[0] << theBand;
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname_opt[index] << " to " << theValue << std::endl;
+ switch( fieldType ){
+ case OFTInteger:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),static_cast<int>(theValue));
+ break;
+ case OFTString:
+ {
+ ostringstream os;
+ os << theValue;
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),os.str().c_str());
+ break;
+ }
+ case OFTReal:
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ else
+ writeCentroidFeature->SetField(fieldname_opt[index].c_str(),theValue);
+ break;
+ case OFTRealList:{
+ int fieldIndex;
+ int nCount;
+ const double *theList;
+ vector<double> vectorList;
+ if(polygon_opt[0]){
+ fieldIndex=writePolygonFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writePolygonFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writePolygonFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ else{
+ fieldIndex=writeCentroidFeature->GetFieldIndex(fieldname_opt[index].c_str());
+ theList=writeCentroidFeature->GetFieldAsDoubleList(fieldIndex,&nCount);
+ vectorList.resize(nCount+1);
+ for(int index=0;index<nCount;++index)
+ vectorList[index]=theList[index];
+ vectorList[nCount]=theValue;
+ writeCentroidFeature->SetField(fieldIndex,vectorList.size(),&(vectorList[0]));
+ }
+ break;
+ }
+ }
+ }
+ break;
+ }
+ case(rule::proportion):{//proportion classes
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyValues);
+ // stat.sum(polyValues);
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ ostringstream fs;
+ fs << class_opt[index];
+ if(polygon_opt[0])
+ writePolygonFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
+ else
+ writeCentroidFeature->SetField(fs.str().c_str(),static_cast<int>(theValue));
+ }
+ break;
+ }
+ case(rule::custom):{//custom
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ stat.normalize_pct(polyValues);
+ assert(polyValues.size()==2);//11:broadleaved, 12:coniferous
+ if(polyValues[0]>=75)//broadleaved
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(11));
+ else if(polyValues[1]>=75)//coniferous
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(12));
+ else if(polyValues[0]>25&&polyValues[1]>25)//mixed
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(13));
+ else{
+ if(verbose_opt[0]){
+ std::cout << "No valid value in polyValues..." << std::endl;
+ for(int index=0;index<polyValues.size();++index){
+ double theValue=polyValues[index];
+ std::cout << theValue << " ";
+ }
+ std::cout << std::endl;
+ }
+ writePolygonFeature->SetField(label_opt[0].c_str(),static_cast<int>(20));
+ }
+ break;
+ }
+ case(rule::minimum):{//minimum of polygon
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for min class
+ int minClass=stat.max(class_opt);
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(polyValues[iclass]>0){
+ if(verbose_opt[0]>1)
+ std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
+ if(class_opt[iclass]<minClass)
+ minClass=class_opt[iclass];
+ }
+ }
+ if(verbose_opt[0]>0)
+ std::cout << "minClass: " << minClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),minClass);
+ break;
+ }
+ case(rule::maximum):{//maximum of polygon
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for max class
+ int maxClass=stat.min(class_opt);
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(polyValues[iclass]>0){
+ if(verbose_opt[0]>1)
+ std::cout << class_opt[iclass] << ": " << polyValues[iclass] << std::endl;
+ if(class_opt[iclass]>maxClass)
+ maxClass=class_opt[iclass];
+ }
+ }
+ if(verbose_opt[0]>0)
+ std::cout << "maxClass: " << maxClass << std::endl;
+ writePolygonFeature->SetField(label_opt[0].c_str(),maxClass);
+ break;
+ }
+ case(rule::maxvote):{//maximum votes in polygon
+ assert(polygon_opt[0]);//not implemented for points
+ if(verbose_opt[0])
+ std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
+ //search for max votes
+ int maxClass=stat.min(class_opt);
+ vector<double>::iterator maxit;
+ maxit=stat.max(polyValues,polyValues.begin(),polyValues.end());
+ int maxIndex=distance(polyValues.begin(),maxit);
+ maxClass=class_opt[maxIndex];
+ }
+ }
+ if(polygon_opt[0]){
+ if(keepFeatures_opt[0])
+ writePolygonFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(verbose_opt[0]>1)
+ std::cout << "creating polygon feature" << std::endl;
if(writeTest){
if(writeTestLayer->CreateFeature( writePolygonFeature ) != OGRERR_NONE ){
std::string errorString="Failed to create polygon feature in shapefile";
@@ -2347,16 +2367,16 @@ int main(int argc, char *argv[])
throw(errorString);
}
}
- OGRFeature::DestroyFeature( writePolygonFeature );
+ OGRFeature::DestroyFeature( writePolygonFeature );
++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
}
- else{
- if(keepFeatures_opt[0])
- writeCentroidFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
- if(verbose_opt[0]>1)
- std::cout << "creating point feature in centroid" << std::endl;
+ else{
+ if(keepFeatures_opt[0])
+ writeCentroidFeature->SetField("origId",static_cast<int>(readFeature->GetFID()));
+ if(verbose_opt[0]>1)
+ std::cout << "creating point feature in centroid" << std::endl;
if(writeTest){
if(writeTestLayer->CreateFeature( writeCentroidFeature ) != OGRERR_NONE ){
std::string errorString="Failed to create point feature in shapefile";
@@ -2369,35 +2389,38 @@ int main(int argc, char *argv[])
throw(errorString);
}
}
- OGRFeature::DestroyFeature( writeCentroidFeature );
+ OGRFeature::DestroyFeature( writeCentroidFeature );
++ntotalvalid;
if(verbose_opt[0])
std::cout << "ntotalvalid: " << ntotalvalid << std::endl;
}
}
}
- else{
- std::string test;
- test=poGeometry->getGeometryName();
- ostringstream oss;
- oss << "geometry " << test << " not supported";
- throw(oss.str());
+ else{
+ std::string test;
+ test=poGeometry->getGeometryName();
+ ostringstream oss;
+ oss << "geometry " << test << " not supported";
+ throw(oss.str());
}
- ++ifeature;
- progress=static_cast<float>(ifeature+1)/nfeature;
- pfnProgress(progress,pszMessage,pProgressArg);
+ ++ifeature;
+ progress=static_cast<float>(ifeature+1)/nfeature;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ catch(std::string e){
+ std::cout << e << std::endl;
+ continue;
}
- catch(std::string e){
- std::cout << e << std::endl;
- continue;
- }
}//end of getNextFeature
- if(rbox_opt[0]>0||cbox_opt[0]>0)
- boxWriter.close();
- ogrWriter.close();
- if(test_opt.size())
- ogrTestWriter.close();
+ // if(rbox_opt[0]>0||cbox_opt[0]>0)
+ // boxWriter.close();
+ progress=1.0;
+ pfnProgress(progress,pszMessage,pProgressArg);
}
+ ogrWriter.close();
+ if(test_opt.size())
+ ogrTestWriter.close();
+ }
progress=1.0;
pfnProgress(progress,pszMessage,pProgressArg);
imgReader.close();
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 9ca1aaf..e24c5f6 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -39,6 +39,7 @@ using namespace std;
int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i","input","input image file");
Optionpk<std::string> output_opt("o", "output", "Output image file");
+ Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
Optionpk<bool> disc_opt("c", "circular", "circular disc kernel for dilation and erosion", false);
// Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation (North=0, East=90, South=180, West=270).");
Optionpk<std::string> method_opt("f", "filter", "filter function (median, var, min, max, sum, mean, dilate, erode, close, open, homog (central pixel must be identical to all other pixels within window), heterog, sobelx (horizontal edge detection), sobely (vertical edge detection), sobelxy (diagonal edge detection NE-SW),sobelyx (diagonal edge detection NW-SE), smooth, density, majority voting (only for classes), smoothnodata (smooth nodata values only) values, threshold local filtering [...]
@@ -75,6 +76,7 @@ int main(int argc,char **argv) {
try{
doProcess=input_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
+ tmpdir_opt.retrieveOption(argc,argv);
disc_opt.retrieveOption(argc,argv);
// angle_opt.retrieveOption(argc,argv);
method_opt.retrieveOption(argc,argv);
@@ -401,7 +403,7 @@ int main(int argc,char **argv) {
exit(1);
}
ostringstream tmps;
- tmps << "/tmp/dilation_" << getpid() << ".tif";
+ tmps << tmpdir_opt[0] << "/dilation_" << getpid() << ".tif";
ImgWriterGdal tmpout;
tmpout.open(tmps.str(),input);
try{
@@ -437,7 +439,7 @@ int main(int argc,char **argv) {
exit(1);
}
ostringstream tmps;
- tmps << "/tmp/erosion_" << getpid() << ".tif";
+ tmps << tmpdir_opt[0] << "/erosion_" << getpid() << ".tif";
ImgWriterGdal tmpout;
tmpout.open(tmps.str(),input);
if(dimZ_opt.size()){
diff --git a/src/apps/pkmosaic.cc b/src/apps/pkmosaic.cc
index d2aa332..39fe661 100644
--- a/src/apps/pkmosaic.cc
+++ b/src/apps/pkmosaic.cc
@@ -46,12 +46,12 @@ int main(int argc, char *argv[])
Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box", 0.0);
Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (0.0: keep original resolution)", 0.0);
Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (0.0: keep original resolution)", 0.0);
- Optionpk<int> band_opt("b", "band", "band index(es) to crop (-1: crop all bands)", -1);
+ Optionpk<int> band_opt("b", "band", "band index(es) to crop (leave empty if all bands must be retained)");
Optionpk<string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
- Optionpk<short> dstnodata_opt("dstnodata", "dstnodata", "nodata value to put in image if out of bounds.", 0);
+ Optionpk<double> dstnodata_opt("dstnodata", "dstnodata", "nodata value to put in image if out of bounds.", 0);
Optionpk<unsigned short> resample_opt("r", "resample", "Resampling method (0: nearest neighbour, 1: bi-linear interpolation).", 0);
Optionpk<string> description_opt("d", "description", "Set image description");
Optionpk<string> crule_opt("cr", "crule", "Composite rule for mosaic (overwrite, maxndvi, maxband, minband, mean, mode (only for byte images), median, sum", "overwrite");
@@ -260,7 +260,7 @@ int main(int argc, char *argv[])
break;
}
}
- if(band_opt[0]>=0){
+ if(band_opt.size()){
nband=band_opt.size();
bands.resize(band_opt.size());
for(int iband=0;iband<band_opt.size();++iband){
@@ -513,7 +513,7 @@ 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[0]<0)?iband:band_opt[iband];
+ int readBand=(band_opt.size()>iband)? band_opt[iband] : iband;
readBuffer[iband].resize(readncol);
try{
imgReader.readData(readBuffer[iband],GDT_Float64,startCol,endCol,readRow,readBand,theResample);
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index 7e51074..21c1da0 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -57,6 +57,26 @@ public:
double getLrx() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lrx);};
double getLry() const {double ulx, uly, lrx,lry;getBoundingBox(ulx,uly,lrx,lry);return(lry);};
// bool getMagicPixel(double& magicX, double& magicY) const {magicX=m_magic_x;magicY=m_magic_y;};
+ void setScale(double theScale, int band=0){
+ /* if(getRasterBand(band)->SetScale(theScale)==CE_Failure){ */
+ if(m_scale.size()!=nrOfBand()){//initialize
+ m_scale.resize(nrOfBand());
+ for(int iband=0;iband<nrOfBand();++iband)
+ m_scale[iband]=1.0;
+ }
+ m_scale[band]=theScale;
+ /* }; */
+ }
+ void setOffset(double theOffset, int band=0){
+ /* if(getRasterBand(band)->SetOffset(theOffset)==CE_Failure){ */
+ if(m_offset.size()!=nrOfBand()){
+ m_offset.resize(nrOfBand());
+ for(int iband=0;iband<nrOfBand();++iband)
+ m_offset[iband]=0.0;
+ }
+ m_offset[band]=theOffset;
+ /* }; */
+ }
int getNoDataValues(std::vector<double>& noDataValues) const;
bool isNoData(double value) const{return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();};
int pushNoDataValue(double noDataValue);
@@ -109,6 +129,8 @@ protected:
/* double m_delta_y; */
bool m_isGeoRef;
std::vector<double> m_noDataValues;
+ std::vector<double> m_scale;
+ std::vector<double> m_offset;
};
// adfGeoTransform[0] /* top left x */
@@ -129,6 +151,10 @@ template<typename T> void ImgReaderGdal::readData(T& value, const GDALDataType&
assert(row<nrOfRow());
assert(row>=0);
poBand->RasterIO(GF_Read,col,row,1,1,&value,1,1,dataType,0,0);
+ if(m_scale.size()>band)
+ value=static_cast<double>(value)*m_scale[band];
+ if(m_offset.size()>band)
+ value=static_cast<double>(value)+m_offset[band];
}
template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const GDALDataType& dataType, int minCol, int maxCol, int row, int band) const
@@ -146,6 +172,16 @@ template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const
if(buffer.size()!=maxCol-minCol+1)
buffer.resize(maxCol-minCol+1);
poBand->RasterIO(GF_Read,minCol,row,buffer.size(),1,&(buffer[0]),buffer.size(),1,dataType,0,0);
+ if(m_scale.size()>band||m_offset.size()>band){
+ double theScale=1;
+ double theOffset=0;
+ if(m_scale.size()>band)
+ theScale=m_scale[band];
+ if(m_offset.size()>band)
+ theOffset=m_offset[band];
+ for(int index=0;index<buffer.size();++index)
+ buffer[index]=theScale*static_cast<double>(buffer[index])+theOffset;
+ }
}
template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, double row, int band, RESAMPLE resample) const
@@ -188,6 +224,12 @@ template<typename T> void ImgReaderGdal::readDataBlock(Vector2d<T>& buffer, cons
template<typename T> void ImgReaderGdal::readDataBlock(std::vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band) const
{
+ double theScale=1;
+ double theOffset=0;
+ if(m_scale.size()>band)
+ theScale=m_scale[band];
+ if(m_offset.size()>band)
+ theOffset=m_offset[band];
//fetch raster band
GDALRasterBand *poBand;
assert(band<nrOfBand()+1);
@@ -203,6 +245,10 @@ template<typename T> void ImgReaderGdal::readDataBlock(std::vector<T>& buffer, c
if(buffer.size()!=(maxRow-minRow+1)*(maxCol-minCol+1))
buffer.resize((maxRow-minRow+1)*(maxCol-minCol+1));
poBand->RasterIO(GF_Read,minCol,minRow,maxCol-minCol+1,maxRow-minRow+1,&(buffer[0]),(maxCol-minCol+1),(maxRow-minRow+1),dataType,0,0);
+ if(m_scale.size()>band||m_offset.size()>band){
+ for(int index=0;index<buffer.size();++index)
+ buffer[index]=theScale*buffer[index]+theOffset;
+ }
}
// template<typename T> void ImgReaderGdal::readDataBlock(vector<T>& buffer, const GDALDataType& dataType , int minCol, int maxCol, int minRow, int maxRow, int band) const
diff --git a/src/imageclasses/ImgWriterOgr.cc b/src/imageclasses/ImgWriterOgr.cc
index 357472c..d5e5de2 100644
--- a/src/imageclasses/ImgWriterOgr.cc
+++ b/src/imageclasses/ImgWriterOgr.cc
@@ -41,8 +41,12 @@ ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderO
{
m_filename=filename;
setCodec(imgReaderOgr.getDriver());
- createLayer(filename,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
- copyFields(imgReaderOgr);
+ int nlayer=imgReaderOgr.getDataSource()->GetLayerCount();
+ for(int ilayer=0;ilayer<nlayer;++ilayer){
+ std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
+ createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
+ copyFields(imgReaderOgr,ilayer);
+ }
}
ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderOgr, bool copyData)
@@ -50,58 +54,66 @@ ImgWriterOgr::ImgWriterOgr(const std::string& filename, ImgReaderOgr& imgReaderO
CPLErrorReset();
m_filename=filename;
setCodec(imgReaderOgr.getDriver());
- createLayer(filename,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
- copyFields(imgReaderOgr);
- if(copyData){
- OGRFeature *poFeature;
- while(true){// (poFeature = imgReaderOgr.getLayer()->GetNextFeature()) != NULL ){
- poFeature = imgReaderOgr.getLayer()->GetNextFeature();
- if( poFeature == NULL )
- break;
- OGRFeature *poDstFeature = NULL;
-
- poDstFeature=createFeature();
-// poDstFeature = OGRFeature::CreateFeature(m_datasource->GetLayer(m_datasource->GetLayerCount()-1)->GetLayerDefn());
- if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
- const char* fmt;
- std::string errorString="Unable to translate feature %d from layer %s.\n";
- fmt=errorString.c_str();
- CPLError( CE_Failure, CPLE_AppDefined,
- fmt,
- poFeature->GetFID(), getLayerName().c_str() );
- // CPLError( CE_Failure, CPLE_AppDefined,
- // "Unable to translate feature %d from layer %s.\n",
- // poFeature->GetFID(), getLayerName().c_str() );
+ int nlayer=imgReaderOgr.getDataSource()->GetLayerCount();
+ for(int ilayer=0;ilayer<nlayer;++ilayer){
+ std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
+ createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
+ copyFields(imgReaderOgr,ilayer);
+ if(copyData){
+ OGRFeature *poFeature;
+ while(true){// (poFeature = imgReaderOgr.getLayer()->GetNextFeature()) != NULL ){
+ poFeature = imgReaderOgr.getLayer(ilayer)->GetNextFeature();
+ if( poFeature == NULL )
+ break;
+ OGRFeature *poDstFeature = NULL;
+
+ poDstFeature=createFeature(ilayer);
+ //todo: check here if SetFrom works (experienced segmentation fault)
+ if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
+ const char* fmt;
+ std::string errorString="Unable to translate feature %d from layer %s.\n";
+ fmt=errorString.c_str();
+ CPLError( CE_Failure, CPLE_AppDefined,
+ fmt,
+ poFeature->GetFID(), getLayerName().c_str() );
+ // CPLError( CE_Failure, CPLE_AppDefined,
+ // "Unable to translate feature %d from layer %s.\n",
+ // poFeature->GetFID(), getLayerName().c_str() );
- OGRFeature::DestroyFeature( poFeature );
- OGRFeature::DestroyFeature( poDstFeature );
- }
- poDstFeature->SetFID( poFeature->GetFID() );
- OGRFeature::DestroyFeature( poFeature );
-
- CPLErrorReset();
- if(createFeature( poDstFeature ) != OGRERR_NONE){
- const char* fmt;
- std::string errorString="Unable to translate feature %d from layer %s.\n";
- fmt=errorString.c_str();
- CPLError( CE_Failure, CPLE_AppDefined,
- fmt,
- poFeature->GetFID(), getLayerName().c_str() );
- OGRFeature::DestroyFeature( poDstFeature );
+ OGRFeature::DestroyFeature( poFeature );
+ OGRFeature::DestroyFeature( poDstFeature );
+ }
+ poDstFeature->SetFID( poFeature->GetFID() );
+ OGRFeature::DestroyFeature( poFeature );
+
+ CPLErrorReset();
+ if(createFeature( poDstFeature,ilayer ) != OGRERR_NONE){
+ const char* fmt;
+ std::string errorString="Unable to translate feature %d from layer %s.\n";
+ fmt=errorString.c_str();
+ CPLError( CE_Failure, CPLE_AppDefined,
+ fmt,
+ poFeature->GetFID(), getLayerName().c_str() );
+ OGRFeature::DestroyFeature( poDstFeature );
+ }
+ OGRFeature::DestroyFeature( poDstFeature );
}
- OGRFeature::DestroyFeature( poDstFeature );
}
}
}
//---------------------------------------------------------------------------
-void ImgWriterOgr::open(const std::string& filename, ImgReaderOgr& imageReader)
+void ImgWriterOgr::open(const std::string& filename, ImgReaderOgr& imgReaderOgr)
{
m_filename=filename;
- setCodec(imageReader.getDriver());
- createLayer(filename,imageReader.getProjection(),imageReader.getGeometryType(),NULL);
- copyFields(imageReader);
+ setCodec(imgReaderOgr.getDriver());
+ int nlayer=imgReaderOgr.getDataSource()->GetLayerCount();
+ for(int ilayer=0;ilayer<nlayer;++ilayer){
+ std::string layername = imgReaderOgr.getLayer(ilayer)->GetName();
+ createLayer(layername,imgReaderOgr.getProjection(),imgReaderOgr.getGeometryType(),NULL);
+ copyFields(imgReaderOgr,ilayer);
+ }
}
void ImgWriterOgr::open(const std::string& filename, const std::string& imageType)
@@ -164,17 +176,23 @@ OGRLayer* ImgWriterOgr::createLayer(const std::string& layername, const std::str
//if points: use wkbPoint
//if no constraints on the types geometry to be written: use wkbUnknown
OGRLayer* poLayer;
+ OGRSpatialReference oSRS;
+
if(theProjection!=""){
- if(theProjection.find("EPSPG:")!=std::string::npos){
- int epsg_code=atoi(theProjection.substr(theProjection.find_first_not_of("EPSG:")).c_str());
- OGRSpatialReference oSRS;
- oSRS.importFromEPSG(epsg_code);
- poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
- }
- else{
- OGRSpatialReference oSRS(theProjection.c_str());
- poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
- }
+ oSRS.SetFromUserInput(theProjection.c_str());
+ poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
+ // if(theProjection.find("EPSPG:")!=std::string::npos){
+ // int epsg_code=atoi(theProjection.substr(theProjection.find_first_not_of("EPSG:")).c_str());
+ // OGRSpatialReference oSRS;
+ // oSRS.importFromEPSG(epsg_code);
+ // poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
+ // }
+ // else{
+ // OGRSpatialReference oSRS(theProjection.c_str());
+ // poLayer=m_datasource->CreateLayer( layername.c_str(), &oSRS, eGType,papszOptions );
+ // }
+ // }
+ // oSRS.importFromProj4(theProjection);
}
else
poLayer=m_datasource->CreateLayer( layername.c_str(), NULL, eGType,papszOptions );
@@ -184,7 +202,7 @@ OGRLayer* ImgWriterOgr::createLayer(const std::string& layername, const std::str
std::string errorstring="Layer creation failed";
throw(errorstring);
}
- OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
+ // OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
return poLayer;
}
@@ -240,7 +258,7 @@ int ImgWriterOgr::getFields(std::vector<OGRFieldDefn*>& fields, int layer) const
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
fields[iField]=poFDefn->GetFieldDefn(iField);
}
- assert(fields.size()==getFieldCount());
+ assert(fields.size()==getFieldCount(layer));
return(fields.size());
}
@@ -262,9 +280,9 @@ void ImgWriterOgr::copyFields(const ImgReaderOgr& imgReaderOgr, int theLayer){
}
}
-void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, const std::string& theId){
+void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, const std::string& theId, int layer){
OGRFeature *poFeature;
- poFeature=createFeature();
+ poFeature=createFeature(layer);
OGRPoint pt;
poFeature->SetField( fieldName.c_str(), theId.c_str());
for(std::map<std::string,double>::const_iterator mit=pointAttributes.begin();mit!=pointAttributes.end();++mit){
@@ -273,16 +291,16 @@ void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,doubl
pt.setX(x);
pt.setY(y);
poFeature->SetGeometry( &pt );
- if(createFeature(poFeature)!=OGRERR_NONE){
+ if(createFeature(poFeature,layer)!=OGRERR_NONE){
std::string errorString="Failed to create feature in shapefile";
throw(errorString);
}
OGRFeature::DestroyFeature( poFeature );
}
-void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, int theId){
+void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, int theId, int layer){
OGRFeature *poFeature;
- poFeature = createFeature();
+ poFeature = createFeature(layer);
OGRPoint pt;
if(pointAttributes.size()+1!=poFeature->GetFieldCount()){
std::ostringstream ess;
@@ -298,7 +316,7 @@ void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,doubl
pt.setX(x);
pt.setY(y);
poFeature->SetGeometry( &pt );
- if(createFeature(poFeature)!=OGRERR_NONE){
+ if(createFeature(poFeature,layer)!=OGRERR_NONE){
std::string errorString="Failed to create feature in shapefile";
throw(errorString);
}
@@ -306,9 +324,9 @@ void ImgWriterOgr::addPoint(double x, double y, const std::map<std::string,doubl
}
//add a line std::string (polygon), caller is responsible to close the line (end point=start point)
-void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId){
+void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer){
OGRFeature *poFeature;
- poFeature = createFeature();
+ poFeature = createFeature(layer);
poFeature->SetStyleString("PEN(c:#FF0000,w:5px)");//see also http://www.gdal.org/ogr/ogr_feature_style.html
poFeature->SetField( fieldName.c_str(), theId);
OGRLineString theLineString;
@@ -319,7 +337,7 @@ void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::stri
std::string errorString="Failed to set line OGRLineString as feature geometry";
throw(errorString);
}
- if(createFeature(poFeature)!=OGRERR_NONE){
+ if(createFeature(poFeature,layer)!=OGRERR_NONE){
std::string errorString="Failed to create feature in shapefile";
throw(errorString);
}
@@ -327,9 +345,9 @@ void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::stri
}
//add a ring (polygon), caller is responsible to close the line (end point=start point)?
-void ImgWriterOgr::addRing(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId){
+void ImgWriterOgr::addRing(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer){
OGRFeature *poFeature;
- poFeature = createFeature();
+ poFeature = createFeature(layer);
poFeature->SetStyleString("PEN(c:#FF0000,w:5px)");//see also http://www.gdal.org/ogr/ogr_feature_style.html
poFeature->SetField( fieldName.c_str(), theId);
// OGRLineString theLineString;
@@ -343,7 +361,7 @@ void ImgWriterOgr::addRing(std::vector<OGRPoint*>& points, const std::string& fi
thePolygon.addRing(&theRing);
// SetSpatialFilter(&thePolygon)
poFeature->SetGeometry( &thePolygon );
- if(createFeature(poFeature)!=OGRERR_NONE){
+ if(createFeature(poFeature,layer)!=OGRERR_NONE){
std::string errorString="Failed to create feature in shapefile";
throw(errorString);
OGRFeature::DestroyFeature( poFeature );
@@ -356,9 +374,9 @@ void ImgWriterOgr::addRing(std::vector<OGRPoint*>& points, const std::string& fi
}
//add a line string (polygon), caller is responsible to close the line (end point=start point)
-void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, const std::string& theId){
+void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, const std::string& theId, int layer){
OGRFeature *poFeature;
- poFeature = createFeature();
+ poFeature = createFeature(layer);
poFeature->SetField( fieldName.c_str(), theId.c_str());
OGRLineString theLineString;
theLineString.setNumPoints(points.size());
@@ -368,25 +386,23 @@ void ImgWriterOgr::addLineString(std::vector<OGRPoint*>& points, const std::stri
std::string errorString="Failed to set line OGRLineString as feature geometry";
throw(errorString);
}
- if(createFeature(poFeature)!=OGRERR_NONE){
+ if(createFeature(poFeature,layer)!=OGRERR_NONE){
std::string errorString="Failed to create feature in shapefile";
throw(errorString);
}
OGRFeature::DestroyFeature( poFeature );
}
-OGRFeature* ImgWriterOgr::createFeature(){
- return(OGRFeature::CreateFeature(m_datasource->GetLayer(m_datasource->GetLayerCount()-1)->GetLayerDefn()));
+OGRFeature* ImgWriterOgr::createFeature(int layer){
+ return(OGRFeature::CreateFeature(m_datasource->GetLayer(layer)->GetLayerDefn()));
}
-OGRErr ImgWriterOgr::createFeature(OGRFeature *theFeature){
- return m_datasource->GetLayer(m_datasource->GetLayerCount()-1)->CreateFeature(theFeature);
+OGRErr ImgWriterOgr::createFeature(OGRFeature *theFeature,int layer){
+ return m_datasource->GetLayer(layer)->CreateFeature(theFeature);
}
int ImgWriterOgr::getFieldCount(int layer) const
{
- if(layer<0)
- layer=m_datasource->GetLayerCount()-1;
assert(m_datasource->GetLayerCount()>layer);
OGRLayer *poLayer;
if((poLayer = m_datasource->GetLayer(layer))==NULL){
@@ -397,9 +413,9 @@ int ImgWriterOgr::getFieldCount(int layer) const
return(poFDefn->GetFieldCount());
}
-int ImgWriterOgr::getFeatureCount() const
+int ImgWriterOgr::getFeatureCount(int layer) const
{
- return(getLayer()->GetFeatureCount());
+ return(getLayer(layer)->GetFeatureCount());
}
int ImgWriterOgr::ascii2ogr(const std::string& filename, const std::string &layername, const std::vector<std::string>& fieldName, const std::vector<OGRFieldType>& fieldType, short colX, short colY, const std::string& theProjection, const OGRwkbGeometryType& eGType, const char fs)
@@ -572,8 +588,6 @@ int ImgWriterOgr::ascii2ogr(const std::string& filename, const std::string &laye
int ImgWriterOgr::addData(const ImgReaderGdal& imgReader, int theLayer, bool verbose)
{
OGRLayer *poLayer;
- if(theLayer<0)
- theLayer=m_datasource->GetLayerCount()-1;//get back layer
assert(m_datasource->GetLayerCount()>theLayer);
if(verbose)
std::cout << "number of layers: " << m_datasource->GetLayerCount() << std::endl;
diff --git a/src/imageclasses/ImgWriterOgr.h b/src/imageclasses/ImgWriterOgr.h
index b48099a..6d54a69 100644
--- a/src/imageclasses/ImgWriterOgr.h
+++ b/src/imageclasses/ImgWriterOgr.h
@@ -46,22 +46,22 @@ public:
int ascii2ogr(const std::string& filename, const std::string &layername, const std::vector<std::string>& fieldName, const std::vector<OGRFieldType>& fieldType, short colX=1, short colY=2, const std::string& theProjection="", const OGRwkbGeometryType& eGType=wkbPoint, const char fs=' ');
OGRLayer* createLayer(const std::string& layername="New layer", const std::string& theProjection="", const OGRwkbGeometryType& eGType=wkbUnknown, char** papszOptions=NULL);
OGRLayer* copyLayer(OGRLayer* poSrcLayer, const std::string& layername, char** papszOptions=NULL);
- void createField(const std::string& fieldname, const OGRFieldType& fieldType, int theLayer=-1);//default: get back layer
+ void createField(const std::string& fieldname, const OGRFieldType& fieldType, int theLayer=0);//default: get back layer
OGRLayer* getLayer(int layer=0) const {return m_datasource->GetLayer(layer);};
std::string getLayerName(int layer=0){return m_datasource->GetLayer(layer)->GetLayerDefn()->GetName();};
int getFields(std::vector<std::string>& fields, int layer=0) const;
int getFields(std::vector<OGRFieldDefn*>& fields, int layer=0) const;
- void copyFields(const ImgReaderOgr& imgReaderOgr, int theLayer=-1);//default: get back layer
- void addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, const std::string& theId);
- void addRing(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId);
- void addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId);
- void addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, const std::string& theId);
- void addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, int theId);
- int addData(const ImgReaderGdal& imgReader, int layer=-1, bool verbose=false);
- OGRFeature* createFeature();
- OGRErr createFeature(OGRFeature* theFeature);
+ void copyFields(const ImgReaderOgr& imgReaderOgr, int theLayer=0);//default: get back layer
+ void addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, const std::string& theId, int layer=0);
+ void addRing(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer=0);
+ void addLineString(std::vector<OGRPoint*>& points, const std::string& fieldName, int theId, int layer=0);
+ void addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, const std::string& theId, int layer=0);
+ void addPoint(double x, double y, const std::map<std::string,double>& pointAttributes, std::string fieldName, int theId, int layer=0);
+ int addData(const ImgReaderGdal& imgReader, int layer=0, bool verbose=false);
+ OGRFeature* createFeature(int layer=0);
+ OGRErr createFeature(OGRFeature* theFeature, int layer=0);
int getFieldCount(int layer=0) const;
- int getFeatureCount() const;
+ int getFeatureCount(int layer=0) const;
OGRDataSource* getDataSource(void) {return m_datasource;};
OGRSFDriver* getDriver(void) const {return m_datasource->GetDriver();};
@@ -69,7 +69,7 @@ protected:
void setCodec(const std::string& imageType);
void setCodec(OGRSFDriver *poDriver);
-// OGRLayer* getLayer(int layer=-1);
+// OGRLayer* getLayer(int layer=0);
std::string m_filename;
OGRDataSource *m_datasource;
--
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