This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository pktools.
commit bd6578cc57291075567f27254cce8884b8717c00
Author: Bas Couwenberg <sebas...@xs4all.nl>
Date: Wed Feb 7 07:59:39 2018 +0100
New upstream version 2.6.7.3+ds
---
CMakeLists.txt | 12 +-
src/apps/pkextractimg.cc | 2 +-
src/apps/pkextractogr.cc | 1037 +++++++++++++++++++++--------------------
test/data/clipped_shp.dbf | Bin 0 -> 791 bytes
test/data/clipped_shp.prj | 1 +
test/data/clipped_shp.qpj | 1 +
test/data/clipped_shp.shp | Bin 0 -> 27100 bytes
test/data/clipped_shp.shx | Bin 0 -> 180 bytes
test/data/test.tif | Bin 0 -> 93155 bytes
test/data/test_clipped.sqlite | Bin 0 -> 45056 bytes
10 files changed, 551 insertions(+), 502 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6d730fe..b334958 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,7 +35,7 @@ set (PROJECT_SOURCE_DIR src)
# The version number.
set (PKTOOLS_VERSION_MAJOR 2)
set (PKTOOLS_VERSION_MINOR 6)
-set (PKTOOLS_VERSION_PATCH 7.2)
+set (PKTOOLS_VERSION_PATCH 7.3)
set (PKTOOLS_VERSION
"${PKTOOLS_VERSION_MAJOR}.${PKTOOLS_VERSION_MINOR}.${PKTOOLS_VERSION_PATCH}")
set (PKTOOLS_PACKAGE_VERSION
"${PKTOOLS_VERSION_MAJOR}.${PKTOOLS_VERSION_MINOR}.${PKTOOLS_VERSION_PATCH}")
set (PKTOOLS_PACKAGE_STRING "PKTOOLS
${PKTOOLS_VERSION_MAJOR}.${PKTOOLS_VERSION_MINOR}.${PKTOOLS_VERSION_PATCH}")
@@ -135,7 +135,6 @@ if(WIN32)
add_definitions(-D_CRT_NONSTDC_NO_WARNING)
add_definitions(-D_SCL_SECURE_NO_WARNINGS)
endif()
-
if(CMAKE_CXX_FLAGS MATCHES "/W[0-4]")
string(REGEX REPLACE "/W[0-4]" "/W4"
CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
@@ -220,7 +219,6 @@ if (BUILD_WITH_LIBLAS)
include_directories(${Boost_INCLUDE_DIRS})
add_definitions("-DBOOST_ALL_NO_LIB")
endif()
-
# include_directories(${BOOST_INCLUDE_DIR})
# if (MSVC)
# set(BOOST_LIBRARIES -LIBPATH:${BOOST_LIB_PATH}
libboost_filesystem-vc100-mt-1_56.lib libboost_system-vc100-mt-1_56.lib)
@@ -348,16 +346,15 @@ if (PKTOOLS_WITH_UTILITIES)
# target_link_libraries (${PKNLOPTUTILITY} ${PKLIBS})
# set_target_properties(${PKNLOPTUTILITY} PROPERTIES FOLDER
utilities)
# endforeach()
- # endif()
+ # endif()
# add_custom_target(utilities DEPENDS ${PKUTILITIES} ${PKLASUTILITIES}
${PKFANNUTILITIES} ${PKNLOPTUTILITIES})
add_custom_target(utilities DEPENDS ${PKUTILITIES} ${PKLASUTILITIES}
${PKFANNUTILITIES})
- set_target_properties(utilities PROPERTIES FOLDER
+ set_target_properties(utilities PROPERTIES FOLDER
phony)
endif()
-
###############################################################################
###############################################################################
@@ -369,7 +366,7 @@ install (FILES "${CMAKE_CURRENT_BINARY_DIR}/pktools-config"
DESTINATION bin PERM
install (FILES "pktools.pc" DESTINATION lib/pkgconfig PERMISSIONS OWNER_READ
OWNER_WRITE GROUP_READ WORLD_READ)
if (PKTOOLS_WITH_UTILITIES)
- install (TARGETS ${PKUTILITIES} DESTINATION bin PERMISSIONS OWNER_READ
OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE)
+ install (TARGETS ${PKUTILITIES} DESTINATION bin PERMISSIONS OWNER_READ
OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE)
if (BUILD_WITH_LIBLAS)
install (TARGETS ${PKLASUTILITIES} DESTINATION bin PERMISSIONS
OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ
WORLD_EXECUTE)
endif(BUILD_WITH_LIBLAS)
@@ -379,7 +376,6 @@ if (PKTOOLS_WITH_UTILITIES)
# if (BUILD_WITH_NLOPT)
# install (TARGETS ${PKNLOPTUTILITIES} DESTINATION bin
PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE
WORLD_READ WORLD_EXECUTE)
# endif(BUILD_WITH_NLOPT)
-
endif(PKTOOLS_WITH_UTILITIES)
###############################################################################
diff --git a/src/apps/pkextractimg.cc b/src/apps/pkextractimg.cc
index 16ce535..d559b70 100644
--- a/src/apps/pkextractimg.cc
+++ b/src/apps/pkextractimg.cc
@@ -90,7 +90,7 @@ using namespace std;
int main(int argc, char *argv[])
{
Optionpk<string> image_opt("i", "input", "Raster input dataset containing
band information");
- Optionpk<string> sample_opt("s", "sample", "OGR vector dataset with features
to be extracted from input data. Output will contain features with input band
information included. Sample image can also be GDAL raster dataset.");
+ Optionpk<string> sample_opt("s", "sample", "Raster dataset with features to
be extracted from input data. Output will contain features with input band
information included.");
Optionpk<string> output_opt("o", "output", "Output sample dataset");
Optionpk<int> class_opt("c", "class", "Class(es) to extract from input
sample image. Leave empty to extract all valid data pixels from sample
dataset");
Optionpk<float> threshold_opt("t", "threshold", "Probability threshold for
selecting samples (randomly). Provide probability in percentage (>0) or
absolute (<0). Use a single threshold per vector sample layer. If using raster
land cover maps as a sample dataset, you can provide a threshold value for each
class (e.g. -t 80 -t 60). Use value 100 to select all pixels for selected
class(es)", 100);
diff --git a/src/apps/pkextractogr.cc b/src/apps/pkextractogr.cc
index 19af55a..3c85766 100644
--- a/src/apps/pkextractogr.cc
+++ b/src/apps/pkextractogr.cc
@@ -36,81 +36,81 @@ along with pktools. If not, see
<http://www.gnu.org/licenses/>.
/******************************************************************************/
/*! \page pkextractogr pkextractogr
- extract pixel values from raster image using a vector dataset sample
-## SYNOPSIS
+ extract pixel values from raster image using a vector dataset sample
+ ## SYNOPSIS
-<code>
+ <code>
Usage: pkextractogr -i input [-s sample | -rand number | -grid size] -o
output
-</code>
+ </code>
-<code>
+ <code>
Options: [-ln layer]* [-c class]* [-t threshold]* [-f format] [-ft
fieldType] [-b band]* [-r rule]*
Advanced options:
[-sband band -eband band]* [-bndnodata band -srcnodata value]* [-tp
threshold] [-buf value [-circ]]
-</code>
-
-\section pkextractogr_description Description
-
-The utility pkextractogr extracts pixel values from an input raster dataset,
based on the locations you provide via a sample file. Alternatively, a random
sample or systematic grid of points can also be extracted. The sample can be a
vector file with points or polygons. In the case of polygons, you can either
extract the values for all raster pixels that are covered by the polygons, or
extract a single value for each polygon such as the centroid, mean, median,
etc. As output, a new copy [...]
-
-A typical usage of pkextractogr is to prepare a training sample for one of the
classifiers implemented in pktools.
-
-\anchor pkextractogr_rules
-
-Overview of the possible extraction rules:
-
-\section pkextractogr_rules Extraction rules:
-
-extraction rule | output features
---------------- | ---------------
-point | extract a single pixel within the polygon or on each point feature
-allpoints | Extract all pixel values covered by the polygon
-centroid | Extract pixel value at the centroid of the polygon
-mean | Extract average of all pixel values within the polygon
-stdev | Extract standard deviation of all pixel values within the polygon
-median | Extract median of all pixel values within the polygon
-min | Extract minimum value of all pixels within the polygon
-max | Extract maximum value of all pixels within the polygon
-sum | Extract sum of the values of all pixels within the polygon
-mode | Extract the mode of classes within the polygon (classes must be set
with the option class)
-proportion | Extract proportion of class(es) within the polygon (classes must
be set with the option class)
-count | Extract count of class(es) within the polygon (classes must be set
with the option class).
-percentile | Extract percentile as defined by option perc (e.g, 95th
percentile of values covered by polygon)
-
-\section pkextractogr_options Options
- - use either `-short` or `--long` options (both `--long=value` and `--long
value` are supported)
- - short option `-h` shows basic options only, long option `--help` shows all
options
-|short|long|type|default|description|
-|-----|----|----|-------|-----------|
- | i | input | std::string | |Raster input dataset
containing band information |
- | s | sample | std::string | |OGR vector dataset
with features to be extracted from input data. Output will contain features
with input band information included |
- | ln | ln | std::string | |Layer name(s) in
sample (leave empty to select all) |
- | rand | random | unsigned int | |Create simple random
sample of points. Provide number of points to generate |
- | grid | grid | double | |Create systematic grid of
points. Provide cell grid size (in projected units, e.g,. m) |
- | o | output | std::string | |Output sample dataset
|
- | c | class | int | |Class(es) to extract from
input sample image. Leave empty to extract all valid data pixels from sample
dataset. Make sure to set classes if rule is set to mode, proportion or count |
- | t | threshold | float | 100 |Probability threshold for
selecting samples (randomly). Provide probability in percentage (>0) or
absolute (<0). Use a single threshold per vector sample layer |
- | perc | perc | double | 95 |Percentile value used for
rule percentile |
- | f | f | std::string | SQLite |Output sample dataset
format |
- | ft | ftype | std::string | Real |Field type (only Real
or Integer) |
- | b | band | int | |Band index(es) to extract (0
based). Leave empty to use all bands |
- | sband | startband | unsigned short | |Start band sequence
number |
- | eband | endband | unsigned short | |End band sequence
number |
- | r | rule | std::string | centroid |Rule how to report
image information per feature (only for vector sample). point (single point or
at centroid if polygon), allpoints (within polygon), centroid, mean, stdev,
median, proportion, count, min, max, mode, sum, percentile. |
- | bndnodata | bndnodata | int | 0 |Band(s) in input image to
check if pixel is valid (used for srcnodata) |
- | srcnodata | srcnodata | double | |Invalid value(s) for
input image |
- | tp | thresholdPolygon | float | |(absolute) threshold for
selecting samples in each polygon |
- | buf | buffer | short | 0 |Buffer for calculating
statistics for point features (in number of pixels) |
- | circ | circular | bool | false |Use a circular disc kernel
buffer (for vector point sample datasets only, use in combination with buffer
option) |
-
-Usage: pkextractogr -i input [-s sample | -rand number | -grid size] -o output
-r rule
-
-
-Examples
-========
-Some examples how to use pkextractogr can be found \ref examples_pkextract
"here"
+ </code>
+
+ \section pkextractogr_description Description
+
+ The utility pkextractogr extracts pixel values from an input raster dataset,
based on the locations you provide via a sample file. Alternatively, a random
sample or systematic grid of points can also be extracted. The sample can be a
vector file with points or polygons. In the case of polygons, you can either
extract the values for all raster pixels that are covered by the polygons, or
extract a single value for each polygon such as the centroid, mean, median,
etc. As output, a new cop [...]
+
+ A typical usage of pkextractogr is to prepare a training sample for one of
the classifiers implemented in pktools.
+
+ \anchor pkextractogr_rules
+
+ Overview of the possible extraction rules:
+
+ \section pkextractogr_rules Extraction rules:
+
+ extraction rule | output features
+ --------------- | ---------------
+ point | extract a single pixel within the polygon or on each point feature
+ allpoints | Extract all pixel values covered by the polygon
+ centroid | Extract pixel value at the centroid of the polygon
+ mean | Extract average of all pixel values within the polygon
+ stdev | Extract standard deviation of all pixel values within the polygon
+ median | Extract median of all pixel values within the polygon
+ min | Extract minimum value of all pixels within the polygon
+ max | Extract maximum value of all pixels within the polygon
+ sum | Extract sum of the values of all pixels within the polygon
+ mode | Extract the mode of classes within the polygon (classes must be set
with the option class)
+ proportion | Extract proportion of class(es) within the polygon (classes
must be set with the option class)
+ count | Extract count of class(es) within the polygon (classes must be set
with the option class).
+ percentile | Extract percentile as defined by option perc (e.g, 95th
percentile of values covered by polygon)
+
+ \section pkextractogr_options Options
+ - use either `-short` or `--long` options (both `--long=value` and `--long
value` are supported)
+ - short option `-h` shows basic options only, long option `--help` shows all
options
+ |short|long|type|default|description|
+ |-----|----|----|-------|-----------|
+ | i | input | std::string | |Raster input dataset
containing band information |
+ | s | sample | std::string | |OGR vector dataset
with features to be extracted from input data. Output will contain features
with input band information included |
+ | ln | ln | std::string | |Layer name(s) in
sample (leave empty to select all) |
+ | rand | random | unsigned int | |Create simple random
sample of points. Provide number of points to generate |
+ | grid | grid | double | |Create systematic grid of
points. Provide cell grid size (in projected units, e.g,. m) |
+ | o | output | std::string | |Output sample dataset
|
+ | c | class | int | |Class(es) to extract from
input sample image. Leave empty to extract all valid data pixels from sample
dataset. Make sure to set classes if rule is set to mode, proportion or count |
+ | t | threshold | float | 100 |Probability threshold for
selecting samples (randomly). Provide probability in percentage (>0) or
absolute (<0). Use a single threshold per vector sample layer |
+ | perc | perc | double | 95 |Percentile value used for
rule percentile |
+ | f | f | std::string | SQLite |Output sample
dataset format |
+ | ft | ftype | std::string | Real |Field type (only Real
or Integer) |
+ | b | band | int | |Band index(es) to extract (0
based). Leave empty to use all bands |
+ | sband | startband | unsigned short | |Start band sequence
number |
+ | eband | endband | unsigned short | |End band sequence
number |
+ | r | rule | std::string | centroid |Rule how to report
image information per feature (only for vector sample). point (single point or
at centroid if polygon), allpoints (within polygon), centroid, mean, stdev,
median, proportion, count, min, max, mode, sum, percentile. |
+ | bndnodata | bndnodata | int | 0 |Band(s) in input image to
check if pixel is valid (used for srcnodata) |
+ | srcnodata | srcnodata | double | |Invalid value(s) for
input image |
+ | tp | thresholdPolygon | float | |(absolute) threshold for
selecting samples in each polygon |
+ | buf | buffer | short | 0 |Buffer for calculating
statistics for point features (in number of pixels) |
+ | circ | circular | bool | false |Use a circular disc kernel
buffer (for vector point sample datasets only, use in combination with buffer
option) |
+
+ Usage: pkextractogr -i input [-s sample | -rand number | -grid size] -o
output -r rule
+
+
+ Examples
+ ========
+ Some examples how to use pkextractogr can be found \ref examples_pkextract
"here"
**/
namespace rule{
@@ -133,8 +133,8 @@ int main(int argc, char *argv[])
Optionpk<string> ogrformat_opt("f", "f", "Output sample dataset
format","SQLite");
Optionpk<string> ftype_opt("ft", "ftype", "Field type (only Real or
Integer)", "Real");
Optionpk<int> band_opt("b", "band", "Band index(es) to extract (0 based).
Leave empty to use all bands");
- Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band
sequence number");
- Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence
number");
+ Optionpk<unsigned short> bstart_opt("sband", "startband", "Start band
sequence number");
+ Optionpk<unsigned short> bend_opt("eband", "endband", "End band sequence
number");
Optionpk<string> rule_opt("r", "rule", "Rule how to report image information
per feature (only for vector sample). point (single point within polygon),
allpoints (all points within polygon), centroid, mean, stdev, median,
proportion, count, min, max, mode, sum, percentile.","centroid");
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);
@@ -189,6 +189,7 @@ int main(int argc, char *argv[])
}
std::map<std::string, rule::RULE_TYPE> ruleMap;
+ std::map<std::string, std::string> fieldMap;
//initialize ruleMap
ruleMap["point"]=rule::point;
ruleMap["centroid"]=rule::centroid;
@@ -205,6 +206,20 @@ int main(int argc, char *argv[])
ruleMap["percentile"]=rule::percentile;
ruleMap["allpoints"]=rule::allpoints;
+ fieldMap["point"]="point";
+ fieldMap["centroid"]="cntrd";
+ fieldMap["mean"]="mean";
+ fieldMap["stdev"]="stdev";
+ fieldMap["median"]="median";
+ fieldMap["proportion"]="prop";
+ fieldMap["count"]="count";
+ fieldMap["min"]="min";
+ fieldMap["max"]="max";
+ fieldMap["custom"]="custom";
+ fieldMap["mode"]="mode";
+ fieldMap["sum"]="sum";
+ fieldMap["percentile"]="perc";
+ fieldMap["allpoints"]="allp";
statfactory::StatFactory stat;
if(srcnodata_opt.size()){
while(srcnodata_opt.size()<bndnodata_opt.size())
@@ -223,11 +238,11 @@ int main(int argc, char *argv[])
// ImgReaderGdal imgReader;
if(image_opt.empty()){
std::cerr << "No image dataset provided (use option -i). Use --help for
help information";
- exit(1);
+ exit(1);
}
if(output_opt.empty()){
std::cerr << "No output dataset provided (use option -o). Use --help for
help information";
- exit(1);
+ exit(1);
}
try{
imgReader.open(image_opt[0],GA_ReadOnly);
@@ -248,17 +263,17 @@ int main(int argc, char *argv[])
try{
if(bstart_opt.size()){
if(bend_opt.size()!=bstart_opt.size()){
- string errorstring="Error: options for start and end band indexes must
be provided as pairs, missing end band";
- throw(errorstring);
+ string errorstring="Error: options for start and end band indexes must
be provided as pairs, missing end band";
+ throw(errorstring);
}
band_opt.clear();
for(int ipair=0;ipair<bstart_opt.size();++ipair){
- if(bend_opt[ipair]<=bstart_opt[ipair]){
- string errorstring="Error: index for end band must be smaller then
start band";
- throw(errorstring);
- }
- for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
- band_opt.push_back(iband);
+ if(bend_opt[ipair]<=bstart_opt[ipair]){
+ string errorstring="Error: index for end band must be smaller then
start band";
+ throw(errorstring);
+ }
+ for(int iband=bstart_opt[ipair];iband<=bend_opt[ipair];++iband)
+ band_opt.push_back(iband);
}
}
}
@@ -266,7 +281,7 @@ int main(int argc, char *argv[])
cerr << error << std::endl;
exit(1);
}
-
+
int nband=(band_opt.size()) ? band_opt.size() : imgReader.nrOfBand();
if(class_opt.size()){
if(nband>1){
@@ -305,7 +320,7 @@ int main(int argc, char *argv[])
exit(1);
break;
}
-
+
const char* pszMessage;
void* pProgressArg=NULL;
GDALProgressFunc pfnProgress=GDALTermProgress;
@@ -348,7 +363,7 @@ int main(int argc, char *argv[])
double theY=uly-static_cast<double>(rand())/(RAND_MAX)*(uly-lry);
pt.setX(theX);
pt.setY(theY);
- pointFeature->SetGeometry( &pt );
+ pointFeature->SetGeometry( &pt );
if(sampleWriterOgr.createFeature(pointFeature) != OGRERR_NONE ){
string errorString="Failed to create feature in vector dataset";
throw(errorString);
@@ -361,7 +376,7 @@ int main(int argc, char *argv[])
}
}
else if(grid_opt.size()){
- //create systematic grid of points
+ //create systematic grid of points
double ulx,uly,lrx,lry;
imgReader.getBoundingBox(ulx,uly,lrx,lry);
if(uly-grid_opt[0]/2<lry&&ulx+grid_opt[0]/2>lrx){
@@ -384,7 +399,7 @@ int main(int argc, char *argv[])
pointFeature=sampleWriterOgr.createFeature();
pt.setX(theX);
pt.setY(theY);
- pointFeature->SetGeometry( &pt );
+ pointFeature->SetGeometry( &pt );
if(sampleWriterOgr.createFeature(pointFeature) != OGRERR_NONE ){
string errorString="Failed to create feature in vector dataset";
throw(errorString);
@@ -429,9 +444,9 @@ int main(int argc, char *argv[])
sampleReaderOgr.getExtent(vectords_ulx,vectords_uly,vectords_lrx,vectords_lry);
bool hasCoverage=((vectords_ulx < imgReader.getLrx())&&(vectords_lrx >
imgReader.getUlx())&&(vectords_lry < imgReader.getUly())&&(vectords_uly >
imgReader.getLry()));
if(!hasCoverage){
- ostringstream ess;
- ess << "No coverage in " << image_opt[0] << " for any layer in " <<
sample_opt[0] << endl;
- throw(ess.str());
+ ostringstream ess;
+ ess << "No coverage in " << image_opt[0] << " for any layer in " <<
sample_opt[0] << endl;
+ throw(ess.str());
}
ogrWriter.open(output_opt[0],ogrformat_opt[0]);
//if class_opt not set, get number of classes from input image for these
rules
@@ -471,7 +486,7 @@ int main(int argc, char *argv[])
cerr << errorString << endl;
exit(1);
}
-
+
//support multiple layers
int nlayerRead=sampleReaderOgr.getDataSource()->GetLayerCount();
int ilayerWrite=0;
@@ -485,11 +500,11 @@ int main(int argc, char *argv[])
string currentLayername=readLayer->GetName();
int layerIndex=ilayer;
if(layer_opt.size()){
- vector<string>::const_iterator
it=find(layer_opt.begin(),layer_opt.end(),currentLayername);
- if(it==layer_opt.end())
- continue;
- else
- layerIndex=it-layer_opt.begin();
+ vector<string>::const_iterator
it=find(layer_opt.begin(),layer_opt.end(),currentLayername);
+ if(it==layer_opt.end())
+ continue;
+ else
+ layerIndex=it-layer_opt.begin();
}
double layer_ulx;
double layer_uly;
@@ -498,7 +513,7 @@ int main(int argc, char *argv[])
sampleReaderOgr.getExtent(layer_ulx,layer_uly,layer_lrx,layer_lry,ilayer);
bool hasCoverage=((layer_ulx < imgReader.getLrx())&&(layer_lrx >
imgReader.getUlx())&&(layer_lry < imgReader.getUly())&&(layer_uly >
imgReader.getLry()));
if(!hasCoverage)
- continue;
+ continue;
//align bounding box to input image
layer_ulx-=fmod(layer_ulx-imgReader.getUlx(),imgReader.getDeltaX());
@@ -536,7 +551,7 @@ int main(int argc, char *argv[])
buffer_opt[0]=1;
}
}
-
+
//extend bounding box with buffer
if(buffer_opt.size()){
layer_uli-=buffer_opt[0];
@@ -570,7 +585,7 @@ int main(int argc, char *argv[])
break;
case OFTReal:
default:
-
imgReader.readDataBlock(readValuesReal[iband],layer_uli,layer_lri,layer_ulj,layer_lrj,theBand);
+
imgReader.readDataBlock(readValuesReal[iband],layer_uli,layer_lri,layer_ulj,layer_lrj,theBand);
break;
}
}
@@ -580,10 +595,10 @@ int main(int argc, char *argv[])
exit(1);
}
-
+
float theThreshold=(threshold_opt.size()==layer_opt.size())?
threshold_opt[layerIndex]: threshold_opt[0];
cout << "processing layer " << currentLayername << endl;
-
+
bool createPolygon=true;
if(find(rule_opt.begin(),rule_opt.end(),"allpoints")!=rule_opt.end())
createPolygon=false;
@@ -591,47 +606,51 @@ int main(int argc, char *argv[])
OGRLayer *writeLayer;
if(createPolygon){
//create polygon
- if(verbose_opt[0])
- std::cout << "create polygons" << std::endl;
- char **papszOptions=NULL;
- writeLayer=ogrWriter.createLayer(readLayer->GetName(),
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);
}
else{
- if(verbose_opt[0])
- std::cout << "create points in layer " << readLayer->GetName() <<
std::endl;
- char **papszOptions=NULL;
+ 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);
+ writeLayer=ogrWriter.createLayer(readLayer->GetName(),
imgReader.getProjection(), wkbPoint, papszOptions);
}
if(verbose_opt[0])
- std::cout << "copy fields from layer " << ilayer << std::flush <<
std::endl;
+ std::cout << "copy fields from layer " << ilayer << std::flush <<
std::endl;
ogrWriter.copyFields(sampleReaderOgr,ilayer,ilayerWrite);
for(int irule=0;irule<rule_opt.size();++irule){
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt.size()) ? band_opt[iband] : iband;
+ ostringstream fs;
+ if(rule_opt.size()>1||nband==1)
+ fs << fieldMap[rule_opt[irule]];
+ if(nband>1)
+ fs << "b" << theBand;
switch(ruleMap[rule_opt[irule]]){
case(rule::proportion):
case(rule::count):{//count for each class
for(int iclass=0;iclass<class_opt.size();++iclass){
- ostringstream fs;
- fs << class_opt[iclass];
- if(nband>1)
- fs << "b" << theBand;
- string fieldname=fs.str();
- ogrWriter.createField(fieldname,fieldType,ilayerWrite);
+ ostringstream fsclass;
+ fsclass << fs.str() << class_opt[iclass];
+ ogrWriter.createField(fsclass.str(),fieldType,ilayerWrite);
}
break;
}
- default:{
- ostringstream fs;
- if(rule_opt.size()>1||nband==1)
- fs << rule_opt[irule];
- if(nband>1)
- fs << "b" << theBand;
- string fieldname=fs.str();
- ogrWriter.createField(fieldname,fieldType,ilayerWrite);
+ case(rule::percentile):{//for each percentile
+ for(int iperc=0;iperc<percentile_opt.size();++iperc){
+ ostringstream fsperc;
+ fsperc << fs.str() << percentile_opt[iperc];
+ ogrWriter.createField(fsperc.str(),fieldType,ilayerWrite);
+ }
+ break;
}
+ default:
+ ogrWriter.createField(fs.str(),fieldType,ilayerWrite);
+ break;
}
}
}
@@ -646,109 +665,109 @@ int main(int argc, char *argv[])
pfnProgress(progress,pszMessage,pProgressArg);
readLayer->ResetReading();
while( (readFeature = readLayer->GetNextFeature()) != NULL ){
- bool validFeature=false;
- if(verbose_opt[0]>2)
- std::cout << "reading feature " << readFeature->GetFID() << std::endl;
- if(theThreshold>0){//percentual value
- double p=static_cast<double>(rand())/(RAND_MAX);
- p*=100.0;
- if(p>theThreshold){
+ bool validFeature=false;
+ if(verbose_opt[0]>2)
+ std::cout << "reading feature " << readFeature->GetFID() <<
std::endl;
+ if(theThreshold>0){//percentual value
+ double p=static_cast<double>(rand())/(RAND_MAX);
+ p*=100.0;
+ if(p>theThreshold){
continue;//do not select for now, go to next feature
- }
- }
- else{//absolute value
- if(threshold_opt.size()==layer_opt.size()){
- if(ntotalvalidLayer>=-theThreshold){
+ }
+ }
+ else{//absolute value
+ if(threshold_opt.size()==layer_opt.size()){
+ if(ntotalvalidLayer>=-theThreshold){
continue;//do not select any more pixels, go to next column
feature
- }
- }
- else{
- if(ntotalvalid>=-theThreshold){
+ }
+ }
+ else{
+ if(ntotalvalid>=-theThreshold){
continue;//do not select any more pixels, go to next column
feature
- }
- }
- }
- if(verbose_opt[0]>2)
- 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 ){
- OGRPoint readPoint = *((OGRPoint *) poGeometry);
-
- double i_centre,j_centre;
+ }
+ }
+ }
+ if(verbose_opt[0]>2)
+ 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 ){
+ OGRPoint readPoint = *((OGRPoint *) poGeometry);
+
+ double i_centre,j_centre;
imgReader.geo2image(readPoint.getX(),readPoint.getY(),i_centre,j_centre);
- //nearest neighbour
- j_centre=static_cast<int>(j_centre);
- i_centre=static_cast<int>(i_centre);
-
- double uli=i_centre-buffer_opt[0];
- double ulj=j_centre-buffer_opt[0];
- double lri=i_centre+buffer_opt[0];
- double lrj=j_centre+buffer_opt[0];
-
- //nearest neighbour
- ulj=static_cast<int>(ulj);
- uli=static_cast<int>(uli);
- lrj=static_cast<int>(lrj);
- lri=static_cast<int>(lri);
-
- //check if j is out of bounds
-
if(static_cast<int>(ulj)<0||static_cast<int>(ulj)>=imgReader.nrOfRow())
- continue;
- //check if j is out of bounds
-
if(static_cast<int>(uli)<0||static_cast<int>(lri)>=imgReader.nrOfCol())
- continue;
-
- OGRPoint ulPoint,urPoint,llPoint,lrPoint;
- double ulx;
+ //nearest neighbour
+ j_centre=static_cast<int>(j_centre);
+ i_centre=static_cast<int>(i_centre);
+
+ double uli=i_centre-buffer_opt[0];
+ double ulj=j_centre-buffer_opt[0];
+ double lri=i_centre+buffer_opt[0];
+ double lrj=j_centre+buffer_opt[0];
+
+ //nearest neighbour
+ ulj=static_cast<int>(ulj);
+ uli=static_cast<int>(uli);
+ lrj=static_cast<int>(lrj);
+ lri=static_cast<int>(lri);
+
+ //check if j is out of bounds
+
if(static_cast<int>(ulj)<0||static_cast<int>(ulj)>=imgReader.nrOfRow())
+ continue;
+ //check if j is out of bounds
+
if(static_cast<int>(uli)<0||static_cast<int>(lri)>=imgReader.nrOfCol())
+ continue;
+
+ OGRPoint ulPoint,urPoint,llPoint,lrPoint;
+ double ulx;
double uly;
- double lrx;
+ double lrx;
double lry;
- OGRPolygon writePolygon;
+ OGRPolygon writePolygon;
OGRPoint writePoint;
- OGRLinearRing writeRing;
- OGRFeature *writePolygonFeature;
-
- int nPointPolygon=0;
- if(createPolygon){
- if(disc_opt[0]){
- double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
- double
radius=buffer_opt[0]*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
- unsigned short nstep = 25;
- for(int i=0;i<nstep;++i){
- OGRPoint aPoint;
-
aPoint.setX(readPoint.getX()+imgReader.getDeltaX()/2.0+radius*cos(2*PI*i/nstep));
-
aPoint.setY(readPoint.getY()-imgReader.getDeltaY()/2.0+radius*sin(2*PI*i/nstep));
- writeRing.addPoint(&aPoint);
- }
- writePolygon.addRing(&writeRing);
- writePolygon.closeRings();
- }
- else{
- double ulx,uly,lrx,lry;
- imgReader.image2geo(uli,ulj,ulx,uly);
- imgReader.image2geo(lri,lrj,lrx,lry);
- ulPoint.setX(ulx-imgReader.getDeltaX()/2.0);
- ulPoint.setY(uly+imgReader.getDeltaY()/2.0);
- lrPoint.setX(lrx+imgReader.getDeltaX()/2.0);
- lrPoint.setY(lry-imgReader.getDeltaY()/2.0);
- urPoint.setX(lrx+imgReader.getDeltaX()/2.0);
- urPoint.setY(uly+imgReader.getDeltaY()/2.0);
- llPoint.setX(ulx-imgReader.getDeltaX()/2.0);
- llPoint.setY(lry-imgReader.getDeltaY()/2.0);
-
- writeRing.addPoint(&ulPoint);
- writeRing.addPoint(&urPoint);
- writeRing.addPoint(&lrPoint);
- writeRing.addPoint(&llPoint);
- writePolygon.addRing(&writeRing);
- writePolygon.closeRings();
- }
+ OGRLinearRing writeRing;
+ OGRFeature *writePolygonFeature;
+
+ int nPointPolygon=0;
+ if(createPolygon){
+ if(disc_opt[0]){
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ double
radius=buffer_opt[0]*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
+ unsigned short nstep = 25;
+ for(int i=0;i<nstep;++i){
+ OGRPoint aPoint;
+
aPoint.setX(readPoint.getX()+imgReader.getDeltaX()/2.0+radius*cos(2*PI*i/nstep));
+
aPoint.setY(readPoint.getY()-imgReader.getDeltaY()/2.0+radius*sin(2*PI*i/nstep));
+ writeRing.addPoint(&aPoint);
+ }
+ writePolygon.addRing(&writeRing);
+ writePolygon.closeRings();
+ }
+ else{
+ double ulx,uly,lrx,lry;
+ imgReader.image2geo(uli,ulj,ulx,uly);
+ imgReader.image2geo(lri,lrj,lrx,lry);
+ ulPoint.setX(ulx-imgReader.getDeltaX()/2.0);
+ ulPoint.setY(uly+imgReader.getDeltaY()/2.0);
+ lrPoint.setX(lrx+imgReader.getDeltaX()/2.0);
+ lrPoint.setY(lry-imgReader.getDeltaY()/2.0);
+ urPoint.setX(lrx+imgReader.getDeltaX()/2.0);
+ urPoint.setY(uly+imgReader.getDeltaY()/2.0);
+ llPoint.setX(ulx-imgReader.getDeltaX()/2.0);
+ llPoint.setY(lry-imgReader.getDeltaY()/2.0);
+
+ writeRing.addPoint(&ulPoint);
+ writeRing.addPoint(&urPoint);
+ writeRing.addPoint(&lrPoint);
+ writeRing.addPoint(&llPoint);
+ writePolygon.addRing(&writeRing);
+ writePolygon.closeRings();
+ }
writePolygonFeature =
OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
if(writePolygonFeature->SetFrom(readFeature)!= OGRERR_NONE)
cerr << "writing feature failed" << std::endl;
@@ -774,33 +793,33 @@ int main(int argc, char *argv[])
valid=valid&&(indexI<imgReader.nrOfCol());
if(valid){
- if(srcnodata_opt.empty())
- validFeature=true;
- else{
- for(int vband=0;vband<bndnodata_opt.size();++vband){
- switch( fieldType ){
- case OFTInteger:{
- int value;
- value=((readValuesInt[vband])[indexJ])[indexI];
- if(value==srcnodata_opt[vband])
- valid=false;
- break;
- }
- case OFTReal:{
- double value;
- value=((readValuesReal[vband])[indexJ])[indexI];
- if(value==srcnodata_opt[vband])
- valid=false;
- break;
- }
- }
- if(!valid)
- continue;
- else
- validFeature=true;
- }
- }
- }
+ if(srcnodata_opt.empty())
+ validFeature=true;
+ else{
+ for(int vband=0;vband<bndnodata_opt.size();++vband){
+ switch( fieldType ){
+ case OFTInteger:{
+ int value;
+ value=((readValuesInt[vband])[indexJ])[indexI];
+ if(value==srcnodata_opt[vband])
+ valid=false;
+ break;
+ }
+ case OFTReal:{
+ double value;
+ value=((readValuesReal[vband])[indexJ])[indexI];
+ if(value==srcnodata_opt[vband])
+ valid=false;
+ break;
+ }
+ }
+ if(!valid)
+ continue;
+ else
+ validFeature=true;
+ }
+ }
+ }
if(valid){
assert(readValuesReal.size()==nband);
for(int iband=0;iband<nband;++iband){
@@ -809,7 +828,7 @@ int main(int argc, char *argv[])
string fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << "centroid";
+ fs << fieldMap["centroid"];
if(nband>1)
fs << "b" << theBand;
fieldname=fs.str();
@@ -873,7 +892,7 @@ int main(int argc, char *argv[])
string fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << "point";
+ fs << fieldMap["point"];
if(nband>1)
fs << "b" << theBand;
fieldname=fs.str();
@@ -897,15 +916,14 @@ int main(int argc, char *argv[])
if(calculateSpatialStatistics||!createPolygon){
Vector2d<double> polyValues;
vector<double> polyClassValues;
-
+
if(class_opt.size()){
polyClassValues.resize(class_opt.size());
//initialize
for(int iclass=0;iclass<class_opt.size();++iclass)
polyClassValues[iclass]=0;
}
- else
- polyValues.resize(nband);
+ polyValues.resize(nband);
OGRPoint thePoint;
for(int j=ulj;j<=lrj;++j){
@@ -987,22 +1005,22 @@ int main(int argc, char *argv[])
}
}
}
- if(valid&&class_opt.size()){
- short value=0;
- switch( fieldType ){
- case OFTInteger:
- value=((readValuesInt[0])[indexJ])[indexI];
- break;
- case OFTReal:
- value=((readValuesReal[0])[indexJ])[indexI];
- break;
- }
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(value==class_opt[iclass])
- polyClassValues[iclass]+=1;
- }
- }
- else if(valid){
+ // if(valid&&class_opt.size()){
+ // short value=0;
+ // switch( fieldType ){
+ // case OFTInteger:
+ // value=((readValuesInt[0])[indexJ])[indexI];
+ // break;
+ // case OFTReal:
+ // value=((readValuesReal[0])[indexJ])[indexI];
+ // break;
+ // }
+ // for(int iclass=0;iclass<class_opt.size();++iclass){
+ // if(value==class_opt[iclass])
+ // polyClassValues[iclass]+=1;
+ // }
+ // }
+ if(valid){
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt.size()) ? band_opt[iband] : iband;
double value=0;
@@ -1014,14 +1032,19 @@ int main(int argc, char *argv[])
value=((readValuesReal[iband])[indexJ])[indexI];
break;
}
-
+ if(!iband&&class_opt.size()){
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(value==class_opt[iclass])
+ polyClassValues[iclass]+=1;
+ }
+ }
if(verbose_opt[0]>1)
std::cout << ": " << value << std::endl;
if(!createPolygon){//write all points within polygon
string fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << "allpoints";
+ fs << fieldMap["allpoints"];
if(nband>1)
fs << "b" << theBand;
fieldname=fs.str();
@@ -1046,7 +1069,7 @@ int main(int argc, char *argv[])
polyValues[iband].push_back(value);
}
}//iband
- }//else (not class_opt.size())
+ }
if(valid&&!createPolygon){
//write feature
if(verbose_opt[0]>1)
@@ -1077,25 +1100,22 @@ int main(int argc, char *argv[])
continue;
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt.size()) ? band_opt[iband] : iband;
- double theValue=0;
- string fieldname;
+ vector<double> theValue;
+ vector<string> fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << rule_opt[irule];
+ fs << fieldMap[rule_opt[irule]];
if(nband>1)
fs << "b" << theBand;
- fieldname=fs.str();
switch(ruleMap[rule_opt[irule]]){
case(rule::proportion):
stat.normalize_pct(polyClassValues);
case(rule::count):{//count for each class
for(int index=0;index<polyClassValues.size();++index){
- theValue=polyClassValues[index];
- ostringstream fs;
- fs << class_opt[index];
- if(nband>1)
- fs << "b" << theBand;
- fieldname=fs.str();
+ theValue.push_back(polyClassValues[index]);
+ ostringstream fsclass;
+ fsclass << fs.str() << class_opt[index];
+ fieldname.push_back(fsclass.str());
}
break;
}
@@ -1111,60 +1131,75 @@ int main(int argc, char *argv[])
maxClass=class_opt[maxIndex];
if(verbose_opt[0]>0)
std::cout << "maxClass: " << maxClass << std::endl;
- theValue=maxClass;
+ theValue.push_back(maxClass);
+ fieldname.push_back(fs.str());
}
case(rule::mean):
- theValue=stat.mean(polyValues[iband]);
+ theValue.push_back(stat.mean(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::median):
- theValue=stat.median(polyValues[iband]);
+ theValue.push_back(stat.median(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::stdev):
- theValue=sqrt(stat.var(polyValues[iband]));
+ theValue.push_back(sqrt(stat.var(polyValues[iband])));
+ fieldname.push_back(fs.str());
break;
- case(rule::percentile):
-
theValue=stat.percentile(polyValues[iband],polyValues[iband].begin(),polyValues[iband].end(),percentile_opt[0]);
+ case(rule::percentile):{
+ for(int iperc=0;iperc<percentile_opt.size();++iperc){
+
theValue.push_back(stat.percentile(polyValues[iband],polyValues[iband].begin(),polyValues[iband].end(),percentile_opt[iperc]));
+ ostringstream fsperc;
+ fsperc << fs.str() << percentile_opt[iperc];
+ fieldname.push_back(fsperc.str());
+ }
break;
+ }
case(rule::sum):
- theValue=stat.sum(polyValues[iband]);
+ theValue.push_back(stat.sum(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::max):
- theValue=stat.mymax(polyValues[iband]);
+ theValue.push_back(stat.mymax(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::min):
- theValue=stat.mymin(polyValues[iband]);
+ theValue.push_back(stat.mymin(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::centroid):
- theValue=polyValues[iband].back();
+ theValue.push_back(polyValues[iband].back());
+ fieldname.push_back(fs.str());
break;
default://not supported
break;
}
-
- switch( fieldType ){
- case OFTInteger:
-
writePolygonFeature->SetField(fieldname.c_str(),static_cast<int>(theValue));
- break;
- case OFTReal:
-
writePolygonFeature->SetField(fieldname.c_str(),theValue);
- break;
- case OFTString:
-
writePolygonFeature->SetField(fieldname.c_str(),type2string<double>(theValue).c_str());
- break;
- default://not supported
- std::string errorString="field type not supported";
- throw(errorString);
- break;
+ for(int ivalue=0;ivalue<theValue.size();++ivalue){
+ switch( fieldType ){
+ case OFTInteger:
+
writePolygonFeature->SetField(fieldname[ivalue].c_str(),static_cast<int>(theValue[ivalue]));
+ break;
+ case OFTReal:
+
writePolygonFeature->SetField(fieldname[ivalue].c_str(),theValue[ivalue]);
+ break;
+ case OFTString:
+
writePolygonFeature->SetField(fieldname[ivalue].c_str(),type2string<double>(theValue[ivalue]).c_str());
+ break;
+ default://not supported
+ std::string errorString="field type not supported";
+ throw(errorString);
+ break;
+ }
}
}
}
}
}
- //test
+ //test
if(createPolygon&&validFeature){
- // if(createPolygon){
+ // if(createPolygon){
//write polygon feature
- //todo: create only in case of valid feature
+ //todo: create only in case of valid feature
if(verbose_opt[0]>1)
std::cout << "creating polygon feature (1)" << std::endl;
if(writeLayer->CreateFeature( writePolygonFeature ) !=
OGRERR_NONE ){
@@ -1176,34 +1211,34 @@ int main(int argc, char *argv[])
++ntotalvalid;
++ntotalvalidLayer;
}
- }
- else{
- OGRPolygon readPolygon;
- OGRMultiPolygon readMultiPolygon;
+ }
+ else{
+ OGRPolygon readPolygon;
+ OGRMultiPolygon readMultiPolygon;
//get envelope
OGREnvelope* psEnvelope=new OGREnvelope();
- if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
- readPolygon = *((OGRPolygon *) poGeometry);
- readPolygon.closeRings();
- readPolygon.getEnvelope(psEnvelope);
- }
- else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon){
- readMultiPolygon = *((OGRMultiPolygon *) poGeometry);
- readMultiPolygon.closeRings();
- readMultiPolygon.getEnvelope(psEnvelope);
- }
- else{
- std::string test;
- test=poGeometry->getGeometryName();
- ostringstream oss;
- oss << "geometry " << test << " not supported";
- throw(oss.str());
- }
-
- double ulx,uly,lrx,lry;
- double uli,ulj,lri,lrj;
+ if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
+ readPolygon = *((OGRPolygon *) poGeometry);
+ readPolygon.closeRings();
+ readPolygon.getEnvelope(psEnvelope);
+ }
+ else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon){
+ readMultiPolygon = *((OGRMultiPolygon *) poGeometry);
+ readMultiPolygon.closeRings();
+ readMultiPolygon.getEnvelope(psEnvelope);
+ }
+ else{
+ std::string test;
+ test=poGeometry->getGeometryName();
+ ostringstream oss;
+ oss << "geometry " << test << " not supported";
+ throw(oss.str());
+ }
+
+ double ulx,uly,lrx,lry;
+ double uli,ulj,lri,lrj;
ulx=psEnvelope->MinX;
uly=psEnvelope->MaxY;
lrx=psEnvelope->MaxX;
@@ -1214,9 +1249,9 @@ int main(int argc, char *argv[])
if(!imgReader.covers(ulx,uly,lrx,lry))
continue;
- OGRFeature *writePolygonFeature;
- int nPointPolygon=0;
- if(createPolygon){
+ OGRFeature *writePolygonFeature;
+ int nPointPolygon=0;
+ if(createPolygon){
writePolygonFeature =
OGRFeature::CreateFeature(writeLayer->GetLayerDefn());
writePolygonFeature->SetGeometry(poGeometry);
//writePolygonFeature and readFeature are both of type wkbPolygon
@@ -1226,16 +1261,16 @@ int main(int argc, char *argv[])
std::cout << "copying new fields write polygon " << std::endl;
if(verbose_opt[0]>1)
std::cout << "write feature has " <<
writePolygonFeature->GetFieldCount() << " fields" << std::endl;
- }
+ }
OGRPoint readPoint;
-
if(find(rule_opt.begin(),rule_opt.end(),"centroid")!=rule_opt.end()){
+
if(find(rule_opt.begin(),rule_opt.end(),"centroid")!=rule_opt.end()){
if(verbose_opt[0]>1)
std::cout << "get centroid" << std::endl;
- if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon)
- readPolygon.Centroid(&readPoint);
- else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon)
- readMultiPolygon.Centroid(&readPoint);
+ if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon)
+ readPolygon.Centroid(&readPoint);
+ else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon)
+ readMultiPolygon.Centroid(&readPoint);
double i,j;
imgReader.geo2image(readPoint.getX(),readPoint.getY(),i,j);
@@ -1247,33 +1282,33 @@ int main(int argc, char *argv[])
valid=valid&&(indexI>=0);
valid=valid&&(indexI<imgReader.nrOfCol());
if(valid){
- if(srcnodata_opt.empty())
- validFeature=true;
- else{
- for(int vband=0;vband<bndnodata_opt.size();++vband){
- switch( fieldType ){
- case OFTInteger:{
- int value;
- value=((readValuesInt[vband])[indexJ])[indexI];
- if(value==srcnodata_opt[vband])
- valid=false;
- break;
- }
- case OFTReal:{
- double value;
- value=((readValuesReal[vband])[indexJ])[indexI];
- if(value==srcnodata_opt[vband])
- valid=false;
- break;
- }
- }
- if(!valid)
- continue;
- else
- validFeature=true;
- }
- }
- }
+ if(srcnodata_opt.empty())
+ validFeature=true;
+ else{
+ for(int vband=0;vband<bndnodata_opt.size();++vband){
+ switch( fieldType ){
+ case OFTInteger:{
+ int value;
+ value=((readValuesInt[vband])[indexJ])[indexI];
+ if(value==srcnodata_opt[vband])
+ valid=false;
+ break;
+ }
+ case OFTReal:{
+ double value;
+ value=((readValuesReal[vband])[indexJ])[indexI];
+ if(value==srcnodata_opt[vband])
+ valid=false;
+ break;
+ }
+ }
+ if(!valid)
+ continue;
+ else
+ validFeature=true;
+ }
+ }
+ }
if(valid){
assert(readValuesReal.size()==nband);
for(int iband=0;iband<nband;++iband){
@@ -1282,7 +1317,7 @@ int main(int argc, char *argv[])
string fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << "centroid";
+ fs << fieldMap["centroid"];
if(nband>1)
fs << "b" << theBand;
fieldname=fs.str();
@@ -1301,17 +1336,17 @@ int main(int argc, char *argv[])
}
}
}
- if(find(rule_opt.begin(),rule_opt.end(),"point")!=rule_opt.end()){
+ if(find(rule_opt.begin(),rule_opt.end(),"point")!=rule_opt.end()){
if(verbose_opt[0]>1)
std::cout << "get point on surface" << std::endl;
- if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
- if(readPolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
- readPolygon.Centroid(&readPoint);
- }
- else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon){
- // if(readMultiPolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
- readMultiPolygon.Centroid(&readPoint);
- }
+ if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
+ if(readPolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
+ readPolygon.Centroid(&readPoint);
+ }
+ else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon){
+ // if(readMultiPolygon.PointOnSurface(&readPoint)!=OGRERR_NONE)
+ readMultiPolygon.Centroid(&readPoint);
+ }
double i,j;
imgReader.geo2image(readPoint.getX(),readPoint.getY(),i,j);
int indexJ=static_cast<int>(j-layer_ulj);
@@ -1352,7 +1387,7 @@ int main(int argc, char *argv[])
string fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << "point";
+ fs << fieldMap["point"];
if(nband>1)
fs << "b" << theBand;
fieldname=fs.str();
@@ -1402,15 +1437,14 @@ int main(int argc, char *argv[])
Vector2d<double> polyValues;
vector<double> polyClassValues;
-
+
if(class_opt.size()){
polyClassValues.resize(class_opt.size());
//initialize
for(int iclass=0;iclass<class_opt.size();++iclass)
polyClassValues[iclass]=0;
}
- else
- polyValues.resize(nband);
+ polyValues.resize(nband);
OGRPoint thePoint;
for(int j=ulj;j<=lrj;++j){
@@ -1429,15 +1463,15 @@ int main(int argc, char *argv[])
thePoint.setX(theX);
thePoint.setY(theY);
//check if point is on surface
- if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
- if(!readPolygon.Contains(&thePoint))
- continue;
- }
- else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon){
- if(!readMultiPolygon.Contains(&thePoint))
- continue;
- }
-
+ if(wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon){
+ if(!readPolygon.Contains(&thePoint))
+ continue;
+ }
+ else if(wkbFlatten(poGeometry->getGeometryType()) ==
wkbMultiPolygon){
+ if(!readMultiPolygon.Contains(&thePoint))
+ continue;
+ }
+
bool valid=true;
if(srcnodata_opt.size()){
for(int vband=0;vband<bndnodata_opt.size();++vband){
@@ -1491,68 +1525,73 @@ int main(int argc, char *argv[])
}
}
}
- if(class_opt.size()){
- short value=0;
+ // if(class_opt.size()){
+ // short value=0;
+ // switch( fieldType ){
+ // case OFTInteger:
+ // value=((readValuesInt[0])[indexJ])[indexI];
+ // break;
+ // case OFTReal:
+ // value=((readValuesReal[0])[indexJ])[indexI];
+ // break;
+ // }
+ // for(int iclass=0;iclass<class_opt.size();++iclass){
+ // if(value==class_opt[iclass])
+ // polyClassValues[iclass]+=1;
+ // }
+ // }
+ // else{
+ for(int iband=0;iband<nband;++iband){
+ int theBand=(band_opt.size()) ? band_opt[iband] : iband;
+ double value=0;
switch( fieldType ){
case OFTInteger:
- value=((readValuesInt[0])[indexJ])[indexI];
+ value=((readValuesInt[iband])[indexJ])[indexI];
break;
case OFTReal:
- value=((readValuesReal[0])[indexJ])[indexI];
+ value=((readValuesReal[iband])[indexJ])[indexI];
break;
}
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(value==class_opt[iclass])
- polyClassValues[iclass]+=1;
+
+ if(!iband&&class_opt.size()){
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(value==class_opt[iclass])
+ polyClassValues[iclass]+=1;
+ }
}
- }
- else{
- for(int iband=0;iband<nband;++iband){
- int theBand=(band_opt.size()) ? band_opt[iband] : iband;
- double value=0;
+ if(verbose_opt[0]>1)
+ std::cout << ": " << value << std::endl;
+ if(!createPolygon){//write all points within polygon
+ string fieldname;
+ ostringstream fs;
+ if(rule_opt.size()>1||nband==1)
+ fs << fieldMap["allpoints"];
+ if(nband>1)
+ fs << "b" << theBand;
+ fieldname=fs.str();
+ int
fieldIndex=writePointFeature->GetFieldIndex(fieldname.c_str());
+ if(fieldIndex<0){
+ cerr << "field " << fieldname << " was not found" <<
endl;
+ exit(1);
+ }
+ if(verbose_opt[0]>1)
+ std::cout << "set field " << fieldname << " to " <<
value << std::endl;
switch( fieldType ){
case OFTInteger:
- value=((readValuesInt[iband])[indexJ])[indexI];
- break;
case OFTReal:
- value=((readValuesReal[iband])[indexJ])[indexI];
+ writePointFeature->SetField(fieldname.c_str(),value);
+ break;
+ default://not supported
+ assert(0);
break;
}
-
- if(verbose_opt[0]>1)
- std::cout << ": " << value << std::endl;
- if(!createPolygon){//write all points within polygon
- string fieldname;
- ostringstream fs;
- if(rule_opt.size()>1||nband==1)
- fs << "allpoints";
- if(nband>1)
- fs << "b" << theBand;
- fieldname=fs.str();
- int
fieldIndex=writePointFeature->GetFieldIndex(fieldname.c_str());
- if(fieldIndex<0){
- cerr << "field " << fieldname << " was not found" <<
endl;
- exit(1);
- }
- if(verbose_opt[0]>1)
- std::cout << "set field " << fieldname << " to " <<
value << std::endl;
- switch( fieldType ){
- case OFTInteger:
- case OFTReal:
- writePointFeature->SetField(fieldname.c_str(),value);
- break;
- default://not supported
- assert(0);
- break;
- }
- }
- else{
- polyValues[iband].push_back(value);
- }
- }//iband
- }//else (not class_opt.size())
+ }
+ else{
+ polyValues[iband].push_back(value);
+ }
+ }//iband
if(!createPolygon){
- //todo: only if valid feature?
+ //todo: only if valid feature?
//write feature
if(verbose_opt[0]>1)
std::cout << "creating point feature" << std::endl;
@@ -1581,25 +1620,22 @@ int main(int argc, char *argv[])
continue;
for(int iband=0;iband<nband;++iband){
int theBand=(band_opt.size()) ? band_opt[iband] : iband;
- double theValue=0;
- string fieldname;
+ vector<double> theValue;
+ vector<string> fieldname;
ostringstream fs;
if(rule_opt.size()>1||nband==1)
- fs << rule_opt[irule];
+ fs << fieldMap[rule_opt[irule]];
if(nband>1)
fs << "b" << theBand;
- fieldname=fs.str();
switch(ruleMap[rule_opt[irule]]){
case(rule::proportion):
stat.normalize_pct(polyClassValues);
case(rule::count):{//count for each class
for(int index=0;index<polyClassValues.size();++index){
- theValue=polyClassValues[index];
- ostringstream fs;
- fs << class_opt[index];
- if(nband>1)
- fs << "b" << theBand;
- fieldname=fs.str();
+ theValue.push_back(polyClassValues[index]);
+ ostringstream fsclass;
+ fsclass << fs.str() << class_opt[index];
+ fieldname.push_back(fsclass.str());
}
break;
}
@@ -1615,58 +1651,73 @@ int main(int argc, char *argv[])
maxClass=class_opt[maxIndex];
if(verbose_opt[0]>0)
std::cout << "maxClass: " << maxClass << std::endl;
- theValue=maxClass;
+ theValue.push_back(maxClass);
+ fieldname.push_back(fs.str());
+ break;
}
case(rule::mean):
- theValue=stat.mean(polyValues[iband]);
+ theValue.push_back(stat.mean(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::median):
- theValue=stat.median(polyValues[iband]);
+ theValue.push_back(stat.median(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::stdev):
- theValue=sqrt(stat.var(polyValues[iband]));
+ theValue.push_back(sqrt(stat.var(polyValues[iband])));
+ fieldname.push_back(fs.str());
break;
- case(rule::percentile):
-
theValue=stat.percentile(polyValues[iband],polyValues[iband].begin(),polyValues[iband].end(),percentile_opt[0]);
+ case(rule::percentile):{
+ for(int iperc=0;iperc<percentile_opt.size();++iperc){
+
theValue.push_back(stat.percentile(polyValues[iband],polyValues[iband].begin(),polyValues[iband].end(),percentile_opt[iperc]));
+ ostringstream fsperc;
+ fsperc << fs.str() << percentile_opt[iperc];
+ fieldname.push_back(fsperc.str());
+ }
break;
+ }
case(rule::sum):
- theValue=stat.sum(polyValues[iband]);
+ theValue.push_back(stat.sum(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::max):
- theValue=stat.mymax(polyValues[iband]);
+ theValue.push_back(stat.mymax(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::min):
- theValue=stat.mymin(polyValues[iband]);
+ theValue.push_back(stat.mymin(polyValues[iband]));
+ fieldname.push_back(fs.str());
break;
case(rule::centroid):
- theValue=polyValues[iband].back();
+ theValue.push_back(polyValues[iband].back());
+ fieldname.push_back(fs.str());
break;
default://not supported
break;
}
-
- switch( fieldType ){
- case OFTInteger:
-
writePolygonFeature->SetField(fieldname.c_str(),static_cast<int>(theValue));
- break;
- case OFTReal:
-
writePolygonFeature->SetField(fieldname.c_str(),theValue);
- break;
- case OFTString:
-
writePolygonFeature->SetField(fieldname.c_str(),type2string<double>(theValue).c_str());
- break;
- default://not supported
- std::string errorString="field type not supported";
- throw(errorString);
- break;
+ for(int ivalue=0;ivalue<theValue.size();++ivalue){
+ switch( fieldType ){
+ case OFTInteger:
+
writePolygonFeature->SetField(fieldname[ivalue].c_str(),static_cast<int>(theValue[ivalue]));
+ break;
+ case OFTReal:
+
writePolygonFeature->SetField(fieldname[ivalue].c_str(),theValue[ivalue]);
+ break;
+ case OFTString:
+
writePolygonFeature->SetField(fieldname[ivalue].c_str(),type2string<double>(theValue[ivalue]).c_str());
+ break;
+ default://not supported
+ std::string errorString="field type not supported";
+ throw(errorString);
+ break;
+ }
}
}
}
}
- }
- //hiero
+ }
if(createPolygon&&validFeature){
- //todo: only create if valid feature?
+ //todo: only create if valid feature?
//write polygon feature
if(verbose_opt[0]>1)
std::cout << "creating polygon feature (2)" << std::endl;
@@ -1679,26 +1730,26 @@ int main(int argc, char *argv[])
++ntotalvalid;
++ntotalvalidLayer;
}
- }
- ++ifeature;
- if(theThreshold>0){
- if(threshold_opt.size()==layer_opt.size())
-
progress=(100.0/theThreshold)*static_cast<float>(ntotalvalidLayer)/nfeatureLayer;
- else
- progress=static_cast<float>(ntotalvalidLayer)/nfeatureLayer;
- }
- else
- progress=static_cast<float>(ifeature+1)/(-theThreshold);
- pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ ++ifeature;
+ if(theThreshold>0){
+ if(threshold_opt.size()==layer_opt.size())
+
progress=(100.0/theThreshold)*static_cast<float>(ntotalvalidLayer)/nfeatureLayer;
+ else
+ progress=static_cast<float>(ntotalvalidLayer)/nfeatureLayer;
+ }
+ else
+ progress=static_cast<float>(ifeature+1)/(-theThreshold);
+ pfnProgress(progress,pszMessage,pProgressArg);
}
- catch(std::string e){
- std::cout << e << std::endl;
- continue;
+ catch(std::string e){
+ std::cout << e << std::endl;
+ continue;
}
- catch(int npoint){
- if(verbose_opt[0])
- std::cout << "number of points read in polygon: " << npoint <<
std::endl;
- continue;
+ catch(int npoint){
+ if(verbose_opt[0])
+ std::cout << "number of points read in polygon: " << npoint <<
std::endl;
+ continue;
}
}
// if(rbox_opt[0]>0||cbox_opt[0]>0)
diff --git a/test/data/clipped_shp.dbf b/test/data/clipped_shp.dbf
new file mode 100644
index 0000000..4223c07
Binary files /dev/null and b/test/data/clipped_shp.dbf differ
diff --git a/test/data/clipped_shp.prj b/test/data/clipped_shp.prj
new file mode 100644
index 0000000..f8ef2c6
--- /dev/null
+++ b/test/data/clipped_shp.prj
@@ -0,0 +1 @@
+PROJCS["IRENET95_Irish_Transverse_Mercator",GEOGCS["GCS_IRENET95",DATUM["D_IRENET95",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",53.5],PARAMETER["central_meridian",-8],PARAMETER["scale_factor",0.99982],PARAMETER["false_easting",600000],PARAMETER["false_northing",750000],UNIT["Meter",1]]
\ No newline at end of file
diff --git a/test/data/clipped_shp.qpj b/test/data/clipped_shp.qpj
new file mode 100644
index 0000000..4f08733
--- /dev/null
+++ b/test/data/clipped_shp.qpj
@@ -0,0 +1 @@
+PROJCS["IRENET95 / Irish Transverse
Mercator",GEOGCS["IRENET95",DATUM["IRENET95",SPHEROID["GRS
1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6173"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4173"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",53.5],PARAMETER["central_meridian",-8],PARAMETER["scale_factor",0.99982],PARAMETER["false_easting
[...]
diff --git a/test/data/clipped_shp.shp b/test/data/clipped_shp.shp
new file mode 100644
index 0000000..c017b5c
Binary files /dev/null and b/test/data/clipped_shp.shp differ
diff --git a/test/data/clipped_shp.shx b/test/data/clipped_shp.shx
new file mode 100644
index 0000000..3a2d10e
Binary files /dev/null and b/test/data/clipped_shp.shx differ
diff --git a/test/data/test.tif b/test/data/test.tif
new file mode 100644
index 0000000..0c3afac
Binary files /dev/null and b/test/data/test.tif differ
diff --git a/test/data/test_clipped.sqlite b/test/data/test_clipped.sqlite
new file mode 100644
index 0000000..247377c
Binary files /dev/null and b/test/data/test_clipped.sqlite differ
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-grass/pktools.git
_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
