[pktools] 03/09: New upstream version 2.6.7.3+ds
Bas Couwenberg
sebastic at debian.org
Wed Feb 7 09:03:10 UTC 2018
This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository pktools.
commit bd6578cc57291075567f27254cce8884b8717c00
Author: Bas Couwenberg <sebastic at 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
More information about the Pkg-grass-devel
mailing list