[pktools] 265/375: removed support for masking in pkexract, introduced srcnodata and bndnodata
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:20 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 15f832e2418dd433fa2144ac890516a8cc2c1ae9
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Sat May 10 23:05:39 2014 +0200
removed support for masking in pkexract, introduced srcnodata and bndnodata
---
ChangeLog | 4 +-
pktools.pc | 2 +-
pktools.pc.in | 2 +-
qt/pkextract_gui/mainwindow.cpp | 15 +-
qt/pkextract_gui/mainwindow.h | 4 -
qt/pkextract_gui/mainwindow.ui | 28 +-
src/apps/pkextract.cc | 730 ++++++++++++++++++++++------------------
src/apps/pkstatogr.cc | 227 +++++++------
8 files changed, 542 insertions(+), 470 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index b798ccf..73e020a 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -293,5 +293,5 @@ version 2.5.2
- pkextract
support for median rule and pointOnSurface
redesign to optimize vector polygon processing
-A
-
+ removed option for masking, introduced srcnodata and bndnodata
+ (more testing is needed)
diff --git a/pktools.pc b/pktools.pc
index 24ad128..fd864e0 100644
--- a/pktools.pc
+++ b/pktools.pc
@@ -7,5 +7,5 @@ Name: pktools
Description: API library for pktools
Requires: gdal gsl
Version: 2.5.2
-Libs: -L${libdir} -lpktools
+Libs: -L${libdir} -lbase -lalgorithms -limageClasses -lfileClasses -llasClasses
Cflags: -I${includedir}/pktools
diff --git a/pktools.pc.in b/pktools.pc.in
index 0731529..55ff1c3 100644
--- a/pktools.pc.in
+++ b/pktools.pc.in
@@ -7,5 +7,5 @@ Name: pktools
Description: API library for pktools
Requires: gdal gsl
Version: @PACKAGE_VERSION@
-Libs: -L${libdir} -lpktools
+Libs: -L${libdir} -lbase -lalgorithms -limageClasses -lfileClasses -llasClasses
Cflags: -I${includedir}/pktools
\ No newline at end of file
diff --git a/qt/pkextract_gui/mainwindow.cpp b/qt/pkextract_gui/mainwindow.cpp
index b1616af..cefc0aa 100644
--- a/qt/pkextract_gui/mainwindow.cpp
+++ b/qt/pkextract_gui/mainwindow.cpp
@@ -30,8 +30,8 @@ void MainWindow::setDefaults()
//tab input/output
ui->input->clear();
ui->sample->clear();
- ui->mask->clear();
- ui->msknodata->setText("0");
+ ui->bndnodata->clear();
+ ui->srcnodata->clear();
ui->polygon->setChecked(false);
ui->f->setCurrentIndex(0);
ui->output->clear();
@@ -56,12 +56,6 @@ void MainWindow::on_actionSample_triggered()
ui->sample->setText(qssample);
}
-void MainWindow::on_actionMask_triggered()
-{
- QString qsmask = QFileDialog::getOpenFileName(this, "Mask");
- ui->mask->setText(qsmask);
-}
-
void MainWindow::on_actionOutput_triggered()
{
QString qsoutput = QFileDialog::getSaveFileName(this,"Output image","","*.*");
@@ -78,11 +72,6 @@ void MainWindow::on_toolButton_sample_clicked()
on_actionSample_triggered();
}
-void MainWindow::on_toolButton_mask_clicked()
-{
- on_actionMask_triggered();
-}
-
void MainWindow::on_toolButton_output_clicked()
{
on_actionOutput_triggered();
diff --git a/qt/pkextract_gui/mainwindow.h b/qt/pkextract_gui/mainwindow.h
index 32e6323..e7e8ce1 100644
--- a/qt/pkextract_gui/mainwindow.h
+++ b/qt/pkextract_gui/mainwindow.h
@@ -20,14 +20,10 @@ private slots:
void on_actionSample_triggered();
- void on_actionMask_triggered();
-
void on_actionOutput_triggered();
void on_toolButton_input_clicked();
- void on_toolButton_mask_clicked();
-
void on_toolButton_output_clicked();
void on_toolButton_sample_clicked();
diff --git a/qt/pkextract_gui/mainwindow.ui b/qt/pkextract_gui/mainwindow.ui
index b6a1372..48fe2f1 100644
--- a/qt/pkextract_gui/mainwindow.ui
+++ b/qt/pkextract_gui/mainwindow.ui
@@ -38,30 +38,23 @@
<item row="1" column="3">
<widget class="QLabel" name="label_14">
<property name="toolTip">
- <string><html><head/><body><p>Pixels in input image where mask image has msknodata will not be classified, but obtain a constant value (=nodata)</p></body></html></string>
+ <string><html><head/><body><p>invalid value(s) for input image</p></body></html></string>
</property>
<property name="text">
- <string>msknodata</string>
+ <string>srcnodata</string>
</property>
</widget>
</item>
<item row="1" column="1">
- <widget class="QLineEdit" name="mask"/>
+ <widget class="QLineEdit" name="bndnodata"/>
</item>
<item row="1" column="0">
<widget class="QLabel" name="label_13">
<property name="toolTip">
- <string><html><head/><body><p>Pixels in input image where mask image has msknodata will not be classified, but obtain a constant value (=nodata)</p></body></html></string>
+ <string><html><head/><body><p>Band(s) in input image to check if pixel is valid (used for srcnodata)</p></body></html></string>
</property>
<property name="text">
- <string>Mask image</string>
- </property>
- </widget>
- </item>
- <item row="1" column="2">
- <widget class="QToolButton" name="toolButton_mask">
- <property name="text">
- <string>...</string>
+ <string>bndnodata</string>
</property>
</widget>
</item>
@@ -104,7 +97,7 @@
<item row="0" column="0">
<widget class="QLabel" name="label_3">
<property name="toolTip">
- <string><html><head/><body><p>Input image (raster or vector) containing band information</p></body></html></string>
+ <string><html><head/><body><p>Raster input dataset containing band information</p></body></html></string>
</property>
<property name="text">
<string>Input data</string>
@@ -117,7 +110,7 @@
<item row="2" column="0">
<widget class="QLabel" name="label_23">
<property name="toolTip">
- <string><html><head/><body><p>OGR vector file with features to be extracted from input data. Output will contain features with input band information included</p></body></html></string>
+ <string><html><head/><body><p>OGR vector file with features to be extracted from input data. Output will contain features with input band information included. Sample image can also be GDAL raster dataset.</p></body></html></string>
</property>
<property name="text">
<string>Sample image</string>
@@ -125,7 +118,7 @@
</widget>
</item>
<item row="1" column="4">
- <widget class="QLineEdit" name="msknodata"/>
+ <widget class="QLineEdit" name="srcnodata"/>
</item>
<item row="3" column="2">
<widget class="QToolButton" name="toolButton_output">
@@ -413,9 +406,8 @@
<tabstops>
<tabstop>input</tabstop>
<tabstop>toolButton_input</tabstop>
- <tabstop>mask</tabstop>
- <tabstop>toolButton_mask</tabstop>
- <tabstop>msknodata</tabstop>
+ <tabstop>bndnodata</tabstop>
+ <tabstop>srcnodata</tabstop>
<tabstop>sample</tabstop>
<tabstop>toolButton_sample</tabstop>
<tabstop>output</tabstop>
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index a6aae69..ffd6a2d 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -42,52 +42,56 @@ using namespace std;
int main(int argc, char *argv[])
{
- Optionpk<string> image_opt("i", "input", "Input image file");
- Optionpk<string> sample_opt("s", "sample", "Input sample vector file or class file (e.g. Corine CLC) if class option is set");
+ Optionpk<string> image_opt("i", "input", "Raster input dataset containing band information");
+ Optionpk<string> sample_opt("s", "sample", "OGR vector file with features to be extracted from input data. Output will contain features with input band information included. Sample image can also be GDAL raster dataset.");
+ Optionpk<string> layer_opt("ln", "ln", "layer name(s) in sample (leave empty to select all)");
Optionpk<string> output_opt("o", "output", "Output sample file (image file)");
Optionpk<int> class_opt("c", "class", "Class(es) to extract from input sample image. Leave empty to extract all valid data pixels from sample file. Make sure to set classes if rule is set to maxvote or proportion");
Optionpk<float> threshold_opt("t", "threshold", "threshold for selecting samples (randomly). Provide probability in percentage (>0) or absolute (<0). Use a single threshold for vector sample files. If using raster land cover maps as a sample file, you can provide a threshold value for each class (e.g. -t 80 -t 60). Use value 100 to select all pixels for selected class(es)", 100);
- Optionpk<float> polythreshold_opt("tp", "thresholdPolygon", "(absolute) threshold for selecting samples in each polygon");
Optionpk<string> ogrformat_opt("f", "f", "Output sample file format","ESRI Shapefile");
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");
Optionpk<bool> polygon_opt("polygon", "polygon", "create OGRPolygon as geometry instead of OGRPoint. Only if sample option is also of polygon type.", false);
- Optionpk<int> band_opt("b", "band", "band index(es) to crop. Use -1 to use all bands)", -1);
- Optionpk<string> rule_opt("r", "rule", "rule how to report image information per feature. point (value at each point or at centroid if polygon), centroid, mean (of polygon), median (of polygon), proportion, minimum (of polygon), maximum (of polygon), maxvote, sum.", "point");
- Optionpk<string> layer_opt("ln", "ln", "layer name(s) in sample (leave empty to select all)");
- Optionpk<string> test_opt("test", "test", "Test sample file (use this option in combination with threshold<100 to create a training (output) and test set");
- Optionpk<string> mask_opt("m", "mask", "Mask image file");
- Optionpk<int> msknodata_opt("msknodata", "msknodata", "Mask value where image is invalid. If a single mask is used, more nodata values can be set. If more masks are used, use one value for each mask.", 1);
+ Optionpk<int> band_opt("b", "band", "band index(es) to extract. Use -1 to use all bands)", -1);
+ Optionpk<string> rule_opt("r", "rule", "rule how to report image information per feature (only for vector sample). point (value at each point or at centroid if polygon), centroid, mean (of polygon), median (of polygon), proportion, minimum (of polygon), maximum (of polygon), maxvote, sum.", "point");
+ Optionpk<double> srcnodata_opt("srcnodata", "srcnodata", "invalid value(s) for input image");
+ Optionpk<int> bndnodata_opt("bndnodata", "bndnodata", "Band(s) in input image to check if pixel is valid (used for srcnodata)", 0);
+ // Optionpk<string> mask_opt("m", "mask", "Mask image file");
+ // Optionpk<int> msknodata_opt("msknodata", "msknodata", "Mask value where image is invalid. If a single mask is used, more nodata values can be set. If more masks are used, use one value for each mask.", 1);
// Optionpk<string> bufferOutput_opt("bu", "bu", "Buffer output shape file");
- Optionpk<string> fieldname_opt("bn", "bname", "Attribute field name of extracted raster band", "B");
+ Optionpk<float> polythreshold_opt("tp", "thresholdPolygon", "(absolute) threshold for selecting samples in each polygon");
+ Optionpk<string> test_opt("test", "test", "Test sample file (use this option in combination with threshold<100 to create a training (output) and test set");
+ Optionpk<string> fieldname_opt("bn", "bname", "For single band input data, this extra attribute name will correspond to the raster values. For multi-band input data, multiple attributes with this prefix will be added (e.g. B0, B1, B2, etc.)", "B");
Optionpk<string> label_opt("cn", "cname", "name of the class label in the output vector file", "label");
Optionpk<short> geo_opt("g", "geo", "use geo coordinates (set to 0 to use image coordinates)", 1);
- Optionpk<short> down_opt("down", "down", "down sampling factor. Can be used to create grid points", 1);
- Optionpk<short> boundary_opt("bo", "boundary", "boundary for selecting the sample", 1);
+ Optionpk<short> down_opt("down", "down", "down sampling factor (for raster sample datasets only). Can be used to create grid points", 1);
+ Optionpk<short> boundary_opt("bo", "boundary", "boundary for selecting the sample (for vector sample datasets only) ", 1);
+ Optionpk<short> disc_opt("circ", "circular", "circular disc kernel boundary (for vector sample datasets only, use in combination with boundary option)", 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<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=image_opt.retrieveOption(argc,argv);
sample_opt.retrieveOption(argc,argv);
+ layer_opt.retrieveOption(argc,argv);
output_opt.retrieveOption(argc,argv);
class_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
- polythreshold_opt.retrieveOption(argc,argv);
ogrformat_opt.retrieveOption(argc,argv);
ftype_opt.retrieveOption(argc,argv);
ltype_opt.retrieveOption(argc,argv);
polygon_opt.retrieveOption(argc,argv);
band_opt.retrieveOption(argc,argv);
rule_opt.retrieveOption(argc,argv);
- layer_opt.retrieveOption(argc,argv);
- test_opt.retrieveOption(argc,argv);
- mask_opt.retrieveOption(argc,argv);
- msknodata_opt.retrieveOption(argc,argv);
+ bndnodata_opt.retrieveOption(argc,argv);
+ srcnodata_opt.retrieveOption(argc,argv);
+ polythreshold_opt.retrieveOption(argc,argv);
+ // mask_opt.retrieveOption(argc,argv);
+ // msknodata_opt.retrieveOption(argc,argv);
// bufferOutput_opt.retrieveOption(argc,argv);
+ test_opt.retrieveOption(argc,argv);
fieldname_opt.retrieveOption(argc,argv);
label_opt.retrieveOption(argc,argv);
geo_opt.retrieveOption(argc,argv);
@@ -120,6 +124,11 @@ int main(int argc, char *argv[])
ruleMap["maxvote"]=rule::maxvote;
ruleMap["sum"]=rule::sum;
+ while(srcnodata_opt.size()<bndnodata_opt.size())
+ srcnodata_opt.push_back(srcnodata_opt[0]);
+ while(bndnodata_opt.size()<srcnodata_opt.size())
+ bndnodata_opt.push_back(bndnodata_opt[0]);
+
if(verbose_opt[0])
std::cout << class_opt << std::endl;
statfactory::StatFactory stat;
@@ -173,26 +182,26 @@ int main(int argc, char *argv[])
if(verbose_opt[0])
std::cout << fieldname_opt << std::endl;
vector<ImgReaderGdal> maskReader;
- if(mask_opt.size()){
- maskReader.resize(mask_opt.size());
- for(int imask=0;imask<mask_opt.size();++imask){
- if(verbose_opt[0]>1)
- std::cout << "opening mask image file " << mask_opt[imask] << std::endl;
- maskReader[imask].open(mask_opt[0]);
- if(imgReader.isGeoRef())
- assert(maskReader[imask].isGeoRef());
- }
- }
+ // if(mask_opt.size()){
+ // maskReader.resize(mask_opt.size());
+ // for(int imask=0;imask<mask_opt.size();++imask){
+ // if(verbose_opt[0]>1)
+ // std::cout << "opening mask image file " << mask_opt[imask] << std::endl;
+ // maskReader[imask].open(mask_opt[0]);
+ // if(imgReader.isGeoRef())
+ // assert(maskReader[imask].isGeoRef());
+ // }
+ // }
- Vector2d<int> maskBuffer;
- if(mask_opt.size()){
- maskBuffer.resize(mask_opt.size());
- for(int imask=0;imask<maskReader.size();++imask)
- maskBuffer[imask].resize(maskReader[imask].nrOfCol());
- }
- vector<double> oldmaskrow(mask_opt.size());
- for(int imask=0;imask<mask_opt.size();++imask)
- oldmaskrow[imask]=-1;
+ // Vector2d<int> maskBuffer;
+ // if(mask_opt.size()){
+ // maskBuffer.resize(mask_opt.size());
+ // for(int imask=0;imask<maskReader.size();++imask)
+ // maskBuffer[imask].resize(maskReader[imask].nrOfCol());
+ // }
+ // vector<double> oldmaskrow(mask_opt.size());
+ // for(int imask=0;imask<mask_opt.size();++imask)
+ // oldmaskrow[imask]=-1;
if(verbose_opt[0]>1)
std::cout << "Number of bands in input image: " << imgReader.nrOfBand() << std::endl;
@@ -343,76 +352,89 @@ int main(int argc, char *argv[])
}
if(static_cast<int>(jimg)<0||static_cast<int>(jimg)>=imgReader.nrOfRow())
continue;
+
+ bool valid=true;
+
if(static_cast<int>(jimg)!=static_cast<int>(oldimgrow)){
assert(imgBuffer.size()==nband);
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt[0]<0)?iband:band_opt[iband];
imgReader.readData(imgBuffer[iband],GDT_Float64,static_cast<int>(jimg),theBand);
assert(imgBuffer[iband].size()==imgReader.nrOfCol());
- }
+ if(bndnodata_opt.size()){
+ vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
+ if(bndit!=bndnodata_opt.end()){
+ vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
+ if(bndit!=bndnodata_opt.end()){
+ if(imgBuffer[iband][static_cast<int>(iimg)]==srcnodata_opt[theBand])
+ valid=false;
+ }
+ }
+ }
+ }
oldimgrow=jimg;
- }
- // 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
- if(geo_opt[0])
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- else{
- colMask=icol;
- rowMask=irow;
- }
- //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;
- }
- }
- 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()){
- if(geo_opt[0])
- maskReader[0].geo2image(x,y,colMask,rowMask);
- else{
- colMask=icol;
- rowMask=irow;
- }
- //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;
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- 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;
- }
- }
- }
- }
+ }
+
+ // for(int imask=0;imask<mask_opt.size();++imask){
+ // double colMask,rowMask;//image coordinates in mask image
+ // if(mask_opt.size()>1){//multiple masks
+ // if(geo_opt[0])
+ // maskReader[imask].geo2image(x,y,colMask,rowMask);
+ // else{
+ // colMask=icol;
+ // rowMask=irow;
+ // }
+ // //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;
+ // }
+ // }
+ // 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()){
+ // if(geo_opt[0])
+ // maskReader[0].geo2image(x,y,colMask,rowMask);
+ // else{
+ // colMask=icol;
+ // rowMask=irow;
+ // }
+ // //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;
+ // if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ // if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ // continue;
+ // 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){
for(int iband=0;iband<imgBuffer.size();++iband){
if(imgBuffer[iband].size()!=imgReader.nrOfCol()){
@@ -587,76 +609,85 @@ int main(int argc, char *argv[])
}
if(static_cast<int>(jimg)<0||static_cast<int>(jimg)>=imgReader.nrOfRow())
continue;
+
+ bool valid=true;
+
if(static_cast<int>(jimg)!=static_cast<int>(oldimgrow)){
assert(imgBuffer.size()==nband);
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt[0]<0)?iband:band_opt[iband];
imgReader.readData(imgBuffer[iband],GDT_Float64,static_cast<int>(jimg),theBand);
assert(imgBuffer[iband].size()==imgReader.nrOfCol());
+ if(bndnodata_opt.size()){
+ vector<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
+ if(bndit!=bndnodata_opt.end()){
+ if(imgBuffer[iband][static_cast<int>(iimg)]==srcnodata_opt[theBand])
+ valid=false;
+ }
+ }
}
oldimgrow=jimg;
}
- 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
- if(geo_opt[0])
- maskReader[imask].geo2image(x,y,colMask,rowMask);
- else{
- colMask=icol;
- rowMask=irow;
- }
- //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;
- }
- }
- 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()){
- if(geo_opt[0])
- maskReader[0].geo2image(x,y,colMask,rowMask);
- else{
- colMask=icol;
- rowMask=irow;
- }
- //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;
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- 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;
- }
- }
- }
- }
+ // for(int imask=0;imask<mask_opt.size();++imask){
+ // double colMask,rowMask;//image coordinates in mask image
+ // if(mask_opt.size()>1){//multiple masks
+ // if(geo_opt[0])
+ // maskReader[imask].geo2image(x,y,colMask,rowMask);
+ // else{
+ // colMask=icol;
+ // rowMask=irow;
+ // }
+ // //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;
+ // }
+ // }
+ // 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()){
+ // if(geo_opt[0])
+ // maskReader[0].geo2image(x,y,colMask,rowMask);
+ // else{
+ // colMask=icol;
+ // rowMask=irow;
+ // }
+ // //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;
+ // if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ // if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ // continue;
+ // 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){
for(int iband=0;iband<imgBuffer.size();++iband){
if(imgBuffer[iband].size()!=imgReader.nrOfCol()){
@@ -819,6 +850,7 @@ int main(int argc, char *argv[])
// assert(fieldnames.size()==ogrWriter.getFieldCount(ilayerWrite));
// map<std::string,double> pointAttributes;
+ //todo: support multiple rules and write attribute for each rule...
if(class_opt.size()){
switch(ruleMap[rule_opt[0]]){
case(rule::proportion):{//proportion for each class
@@ -900,72 +932,70 @@ int main(int argc, char *argv[])
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());
- }
+
+ // 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;
- else
- validFeature=true;
+ // 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;
+ // }
+ // }
+ // }
+ // }
double value;
double i_centre,j_centre;
@@ -1063,6 +1093,20 @@ int main(int argc, char *argv[])
if(verbose_opt[0]>1)
std::cout << "write feature has " << writeFeature->GetFieldCount() << " fields" << std::endl;
+ // //hiero
+ // for(int vband=0;vband<bndnodata_opt.size();++vband){
+ // value=((readValues[bndnodata_opt[vband]])[j-ulj])[i-uli];
+ // if(value==srcnodata_opt[vband]){
+ // valid=false;
+ // break;
+ // }
+ // }
+
+ // if(!valid)
+ // continue;
+ // else
+ // validFeature=true;
+
vector<double> windowBuffer;
for(int windowJ=-theDim/2;windowJ<(theDim+1)/2;++windowJ){
for(int windowI=-theDim/2;windowI<(theDim+1)/2;++windowI){
@@ -1081,6 +1125,15 @@ int main(int argc, char *argv[])
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(bndnodata_opt.size()){
+ Optionpk<int>::const_iterator bndit=find(bndnodata_opt.begin(),bndnodata_opt.end(),theBand);
+ if(bndit!=bndnodata_opt.end()){
+ if(value==srcnodata_opt[theBand])
+ valid=false;
+ }
+ }
+
if(verbose_opt[0]>1)
std::cout << ": " << value << std::endl;
ostringstream fs;
@@ -1258,7 +1311,7 @@ int main(int argc, char *argv[])
assert(lrj<imgReader.nrOfRow());
imgReader.readDataBlock(readValues[iband],GDT_Float64,uli,lri,ulj,lrj,theBand);
}
- //todo: readDataBlock for maskReader...
+
OGRPoint thePoint;
for(int j=ulj;j<=lrj;++j){
for(int i=uli;i<=lri;++i){
@@ -1276,60 +1329,70 @@ int main(int argc, char *argv[])
if(ruleMap[rule_opt[0]]!=rule::centroid&&!readPolygon.Contains(&thePoint))
continue;
+
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());
- }
- }
- 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;
-
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- 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;
- }
- }
+
+ for(int vband=0;vband<bndnodata_opt.size();++vband){
+ double value=((readValues[bndnodata_opt[vband]])[j-ulj])[i-uli];
+ if(value==srcnodata_opt[vband]){
+ valid=false;
+ break;
}
}
+
+ // 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());
+ // }
+ // }
+ // 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;
+
+ // if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ // if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ // continue;
+ // 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;
else
@@ -1828,60 +1891,71 @@ int main(int argc, char *argv[])
if(ruleMap[rule_opt[0]]!=rule::centroid&&!readPolygon.Contains(&thePoint))
continue;
+
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;
+
+ for(int vband=0;vband<bndnodata_opt.size();++vband){
+ double value=((readValues[bndnodata_opt[vband]])[j-ulj])[i-uli];
+ if(value==srcnodata_opt[vband]){
+ valid=false;
+ break;
+ }
+ }
+
+ // 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());
- }
- }
- 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;
+ // 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());
+ // }
+ // }
+ // 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;
- if(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
- if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
- continue;
- 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(static_cast<int>(rowMask)!=static_cast<int>(oldmaskrow[0])){
+ // if(static_cast<int>(rowMask)<0||static_cast<int>(rowMask)>=maskReader[0].nrOfRow())
+ // continue;
+ // 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;
else
diff --git a/src/apps/pkstatogr.cc b/src/apps/pkstatogr.cc
index 5e8b43e..11c34eb 100644
--- a/src/apps/pkstatogr.cc
+++ b/src/apps/pkstatogr.cc
@@ -25,52 +25,56 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include "imageclasses/ImgReaderOgr.h"
#include "algorithms/StatFactory.h"
+using namespace std;
+
int main(int argc, char *argv[])
{
- Optionpk<std::string> input_opt("i", "input", "Input shape file", "");
- Optionpk<std::string> fieldname_opt("n", "fname", "fields on which to calculate statistics", "");
+ Optionpk<string> input_opt("i", "input", "Input OGR vector file", "");
+ Optionpk<string> layer_opt("ln", "lname", "layer name(s) in sample (leave empty to select all)");
+ Optionpk<string> fieldname_opt("n", "fname", "fields on which to calculate statistics", "");
+ Optionpk<double> nodata_opt("nodata","nodata","set nodata value(s)");
+ Optionpk<double> src_min_opt("src_min","src_min","set minimum value for histogram");
+ Optionpk<double> src_max_opt("src_max","src_max","set maximum value for histogram");
+ Optionpk<bool> size_opt("s","size","sample size (number of points)",false);
Optionpk<bool> minmax_opt("mm","minmax","calculate minimum and maximum value",false);
Optionpk<bool> min_opt("min","min","calculate minimum value",0);
Optionpk<bool> max_opt("max","max","calculate maximum value",0);
- Optionpk<double> src_min_opt("src_min","src_min","set minimum value for histogram");
- Optionpk<double> src_max_opt("src_max","src_max","set maximum value for histogram");
- Optionpk<double> nodata_opt("nodata","nodata","set nodata value(s)");
+ Optionpk<bool> mean_opt("mean","mean","calculate mean value",false);
+ Optionpk<bool> median_opt("median","median","calculate median value",false);
+ Optionpk<bool> stdev_opt("stdev","stdev","calculate standard deviation",false);
Optionpk<bool> histogram_opt("hist","hist","calculate histogram",false);
Optionpk<unsigned int> nbin_opt("nbin", "nbin", "number of bins");
Optionpk<bool> relative_opt("rel","relative","use percentiles for histogram to calculate histogram",false);
Optionpk<double> kde_opt("kde","kde","bandwith of kernel density when producing histogram, use 0 for practical estimation based on Silverman's rule of thumb. Leave empty if no kernel density is required");
- Optionpk<bool> mean_opt("mean","mean","calculate mean value",false);
- Optionpk<bool> median_opt("median","median","calculate median value",false);
- Optionpk<bool> stdev_opt("stdev","stdev","calculate standard deviation",false);
- Optionpk<bool> size_opt("s","size","sample size (number of points)",false);
Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
bool doProcess;//stop process when program was invoked with help option (-h --help)
try{
doProcess=input_opt.retrieveOption(argc,argv);
+ layer_opt.retrieveOption(argc,argv);
fieldname_opt.retrieveOption(argc,argv);
+ nodata_opt.retrieveOption(argc,argv);
+ src_min_opt.retrieveOption(argc,argv);
+ src_max_opt.retrieveOption(argc,argv);
+ size_opt.retrieveOption(argc,argv);
minmax_opt.retrieveOption(argc,argv);
min_opt.retrieveOption(argc,argv);
max_opt.retrieveOption(argc,argv);
- src_min_opt.retrieveOption(argc,argv);
- src_max_opt.retrieveOption(argc,argv);
- nodata_opt.retrieveOption(argc,argv);
+ mean_opt.retrieveOption(argc,argv);
+ median_opt.retrieveOption(argc,argv);
+ stdev_opt.retrieveOption(argc,argv);
histogram_opt.retrieveOption(argc,argv);
nbin_opt.retrieveOption(argc,argv);
relative_opt.retrieveOption(argc,argv);
kde_opt.retrieveOption(argc,argv);
- mean_opt.retrieveOption(argc,argv);
- median_opt.retrieveOption(argc,argv);
- stdev_opt.retrieveOption(argc,argv);
- size_opt.retrieveOption(argc,argv);
verbose_opt.retrieveOption(argc,argv);
}
- catch(std::string predefinedString){
- std::cout << predefinedString << std::endl;
+ catch(string predefinedString){
+ cout << predefinedString << endl;
exit(0);
}
if(!doProcess){
- std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
+ cout << "short option -h shows basic options only, use long option --help to show all options" << endl;
exit(0);//help was invoked, stop processing
}
@@ -78,105 +82,122 @@ int main(int argc, char *argv[])
try{
imgReader.open(input_opt[0]);
}
- catch(std::string errorstring){
- std::cerr << errorstring << std::endl;
+ catch(string errorstring){
+ cerr << errorstring << endl;
}
ImgReaderOgr inputReader(input_opt[0]);
- std::vector<double> theData;
+ vector<double> theData;
statfactory::StatFactory stat;
//todo: implement ALL
stat.setNoDataValues(nodata_opt);
- for(int ifield=0;ifield<fieldname_opt.size();++ifield){
+
+ //support multiple layers
+ int nlayerRead=inputReader.getDataSource()->GetLayerCount();
+ if(verbose_opt[0])
+ cout << "number of layers: " << nlayerRead << endl;
+
+ for(int ilayer=0;ilayer<nlayerRead;++ilayer){
+ OGRLayer *readLayer=inputReader.getLayer(ilayer);
+ string currentLayername=readLayer->GetName();
+ if(layer_opt.size())
+ if(find(layer_opt.begin(),layer_opt.end(),currentLayername)==layer_opt.end())
+ continue;
if(verbose_opt[0])
- std::cout << "field: " << ifield << std::endl;
- theData.clear();
- inputReader.readData(theData,OFTReal,fieldname_opt[ifield],0,verbose_opt[0]);
- std::vector<double> binData;
- double minValue=0;
- double maxValue=0;
- stat.minmax(theData,theData.begin(),theData.end(),minValue,maxValue);
- if(src_min_opt.size())
- minValue=src_min_opt[0];
- if(src_max_opt.size())
- maxValue=src_max_opt[0];
- unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
+ cout << "processing layer " << currentLayername << endl;
+ cout << " --lname " << currentLayername;
+
+ for(int ifield=0;ifield<fieldname_opt.size();++ifield){
+ if(verbose_opt[0])
+ cout << "field: " << ifield << endl;
+ theData.clear();
+ inputReader.readData(theData,OFTReal,fieldname_opt[ifield],ilayer,verbose_opt[0]);
+ vector<double> binData;
+ double minValue=0;
+ double maxValue=0;
+ stat.minmax(theData,theData.begin(),theData.end(),minValue,maxValue);
+ if(src_min_opt.size())
+ minValue=src_min_opt[0];
+ if(src_max_opt.size())
+ maxValue=src_max_opt[0];
+ unsigned int nbin=(nbin_opt.size())? nbin_opt[0]:0;
- if(histogram_opt[0]){
- double sigma=0;
- if(kde_opt.size()){
- if(kde_opt[0]>0)
- sigma=kde_opt[0];
- else
- sigma=1.06*sqrt(stat.var(theData))*pow(theData.size(),-0.2);
+ if(histogram_opt[0]){
+ double sigma=0;
+ if(kde_opt.size()){
+ if(kde_opt[0]>0)
+ sigma=kde_opt[0];
+ else
+ sigma=1.06*sqrt(stat.var(theData))*pow(theData.size(),-0.2);
+ }
+ if(nbin<1)
+ nbin=(maxValue-minValue+1);
+ try{
+ stat.distribution(theData,theData.begin(),theData.end(),binData,nbin,minValue,maxValue,sigma);
+ }
+ catch(string theError){
+ cerr << "Warning: all identical values in data" << endl;
+ exit(1);
+ }
}
- if(nbin<1)
- nbin=(maxValue-minValue+1);
+ // int nbin=(nbin_opt[0]>1)? nbin_opt[0] : 2;
+ cout << " --fname " << fieldname_opt[ifield];
try{
- stat.distribution(theData,theData.begin(),theData.end(),binData,nbin,minValue,maxValue,sigma);
- }
- catch(std::string theError){
- std::cerr << "Warning: all identical values in data" << std::endl;
- exit(1);
- }
- }
- // int nbin=(nbin_opt[0]>1)? nbin_opt[0] : 2;
- std::cout << " --fname " << fieldname_opt[ifield];
- try{
- double theMean=0;
- double theVar=0;
- stat.meanVar(theData,theMean,theVar);
- if(mean_opt[0])
- std::cout << " --mean " << theMean;
- if(stdev_opt[0])
- std::cout << " --stdev " << sqrt(theVar);
- if(minmax_opt[0]||min_opt[0]||max_opt[0]){
- if(minmax_opt[0])
- std::cout << " --min " << minValue << " --max " << maxValue << " ";
- else{
- if(min_opt[0])
- std::cout << " --min " << minValue << " ";
- if(max_opt[0])
- std::cout << " --max " << maxValue << " ";
+ double theMean=0;
+ double theVar=0;
+ stat.meanVar(theData,theMean,theVar);
+ if(mean_opt[0])
+ cout << " --mean " << theMean;
+ if(stdev_opt[0])
+ cout << " --stdev " << sqrt(theVar);
+ if(minmax_opt[0]||min_opt[0]||max_opt[0]){
+ if(minmax_opt[0])
+ cout << " --min " << minValue << " --max " << maxValue << " ";
+ else{
+ if(min_opt[0])
+ cout << " --min " << minValue << " ";
+ if(max_opt[0])
+ cout << " --max " << maxValue << " ";
+ }
+ }
+ if(median_opt[0])
+ cout << " -median " << stat.median(theData);
+ if(size_opt[0])
+ cout << " -size " << theData.size();
+ cout << endl;
+ if(histogram_opt[0]){
+ for(int ibin=0;ibin<nbin;++ibin){
+ double binValue=0;
+ if(nbin==maxValue-minValue+1)
+ binValue=minValue+ibin;
+ else
+ binValue=minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin;
+ cout << binValue << " ";
+ if(relative_opt[0])
+ cout << 100.0*static_cast<double>(binData[ibin])/theData.size() << endl;
+ else
+ cout << binData[ibin] << endl;
+ }
}
}
- if(median_opt[0])
- std::cout << " -median " << stat.median(theData);
- if(size_opt[0])
- std::cout << " -size " << theData.size();
- std::cout << std::endl;
- if(histogram_opt[0]){
- for(int ibin=0;ibin<nbin;++ibin){
- double binValue=0;
- if(nbin==maxValue-minValue+1)
- binValue=minValue+ibin;
- else
- binValue=minValue+static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin;
- std::cout << binValue << " ";
- if(relative_opt[0])
- std::cout << 100.0*static_cast<double>(binData[ibin])/theData.size() << std::endl;
- else
- std::cout << binData[ibin] << std::endl;
- }
+ catch(string theError){
+ if(mean_opt[0])
+ cout << " --mean " << theData.back();
+ if(stdev_opt[0])
+ cout << " --stdev " << "0";
+ if(min_opt[0])
+ cout << " -min " << theData.back();
+ if(max_opt[0])
+ cout << " -max " << theData.back();
+ if(median_opt[0])
+ cout << " -median " << theData.back();
+ if(size_opt[0])
+ cout << " -size " << theData.size();
+ cout << endl;
+ cerr << "Warning: all identical values in data" << endl;
}
}
- catch(std::string theError){
- if(mean_opt[0])
- std::cout << " --mean " << theData.back();
- if(stdev_opt[0])
- std::cout << " --stdev " << "0";
- if(min_opt[0])
- std::cout << " -min " << theData.back();
- if(max_opt[0])
- std::cout << " -max " << theData.back();
- if(median_opt[0])
- std::cout << " -median " << theData.back();
- if(size_opt[0])
- std::cout << " -size " << theData.size();
- std::cout << std::endl;
- std::cerr << "Warning: all identical values in data" << std::endl;
- }
}
imgReader.close();
}
--
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