[pktools] 333/375: buffer for points in pkextract
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:27 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit 8fc7db8b8052c264777a097ca6ddd46b282616bd
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Mon Sep 29 18:24:40 2014 +0200
buffer for points in pkextract
---
doc/examples_pkextract.dox | 29 ++++++++++++----
doc/examples_pkfilter.dox | 11 ++++++
src/apps/pkextract.cc | 87 ++++++++++++++++++++++++++++------------------
3 files changed, 87 insertions(+), 40 deletions(-)
diff --git a/doc/examples_pkextract.dox b/doc/examples_pkextract.dox
index d3a9d8d..dfdd665 100644
--- a/doc/examples_pkextract.dox
+++ b/doc/examples_pkextract.dox
@@ -2,36 +2,51 @@
\code
pkextract -i input.tif -s points.sqlite -o extracted.sqlite
\endcode
-extract the points read in points.sqlite from input.tif. Create a new point vector file extracted.sqlite, where each point will contain an attribute for the individual input bands in input.tif.
+Extract all points for all layers read in points.sqlite from input.tif. Create a new point vector dataset named extracted.sqlite, where each point will contain an attribute for the individual input bands in input.tif. Notice that the default vector format is Spatialite (.sqlite).
+
+\code
+pkextract -i input.tif -s points.sqlite -ln valid -o extracted.sqlite
+\endcode
+Same example as above, but only extract the points for the layer in points.sqlite named "valid"
\code
pkextract -i input.tif -s points.shp -f "ESRI Shapefile" -o extracted.shp
\endcode
-Same example as above, but vector format is ESRI Shapefile
+Extract points and write output in ESRI Shapefile format
+
+\code
+pkextract -i input.tif -s points.sqlite -o extracted.sqlite -r stdev -buf 3 -polygon
+\endcode
+Extract the standard deviation for each input band in a 3 by 3 window, centered around the points in the sample vector dataset points.sqlite. The output vector dataset will contain polygon features defined by the buffered points (3x3 window). Use the option -circ to define a circular buffer.
\code
pkextract -i input.tif -s polygons.sqlite -o training.sqlite -r point
\endcode
-extract all pixels from input.tif covered by the polygons in locations.sqlite. Each polygon can thus result in multiple point features with attributes for each input band. Write the extracted points to a point vector file training.sqlite.
+Extract all pixels from input.tif covered by the polygons in locations.sqlite. Each polygon can thus result in multiple point features with attributes for each input band. Write the extracted points to a point vector dataset training.sqlite.
\code
pkextract -i input.tif -b 0 -s polygons.sqlite -r centroid -o extracted.sqlite -polygon
\endcode
-extract the first band from input.tif at the centroids of the polygons in vector filepolygons.sqlite. Assign the extracted point value to a new attribute of the polygon and write to the vector file extracted.sqlite.
+Extract the first band from input.tif at the centroids of the polygons in vector dataset polygons.sqlite. Assign the extracted point value to a new attribute of the polygon and write to the vector dataset extracted.sqlite.
\code
pkextract -i input.tif -b 1 -s polygons.sqlite -r mean -o extracted.sqlite -polygon
\endcode
-extract the mean values for the second band in input.tif covered by each polygon in polygons.sqlite. The mean values are written to a copy of the polygons in output vector file extracted.sqlite
+Extract the mean values for the second band in input.tif covered by each polygon in polygons.sqlite. The mean values are written to a copy of the polygons in output vector dataset extracted.sqlite
+
+\code
+pkextract -i landcover.tif -s polygons.sqlite -r maxvote -o majority.sqlite -polygon -c 1 -c 2 -c 3 -c 4 -c 5
+\endcode
+Extract the majority class in each polygon for the input land cover map. The land cover map contains five valid classes, labeled 1-5. Other class values (e.g., labeled as 0) are not taken into account in the voting.
\code
pkextract -i input.tif -s sample.tif -o extracted.sqlite -t 10 -c 1 -c 2 -c 3
\endcode
-Typical use where pixels are extracted based on a land cover map (sample.tif). Extract all bands for a random sample of 10 percent of the pixels in the land cover map sample.tif where the land cover classes are either 1,2 or 3 (class values). Write output to the point vector file extracted.sqlite.
+Typical use where pixels are extracted based on a land cover map (sample.tif). Extract all bands for a random sample of 10 percent of the pixels in the land cover map sample.tif where the land cover classes are either 1,2 or 3 (class values). Write output to the point vector dataset extracted.sqlite.
\code
pkextract -i input.tif -s sample.tif -o extracted.sqlite -t -5000 -c 1
\endcode
-extract all bands for the first 5000 pixels encountered in sample.tif where pixels have a value equal to 1. Write output to point vector file extracted.sqlite.
+Extract all bands for the first 5000 pixels encountered in sample.tif where pixels have a value equal to 1. Write output to point vector dataset extracted.sqlite.
diff --git a/doc/examples_pkfilter.dox b/doc/examples_pkfilter.dox
index fdb9dcb..7515b23 100644
--- a/doc/examples_pkfilter.dox
+++ b/doc/examples_pkfilter.dox
@@ -8,3 +8,14 @@ filter input.tif with morphological dilation filter. Use a circular kernel (inst
pkfilter -i input.tif -o filter.tif -dx 3 -dy 3 -class 255 -f dilate -c
\endcode
Similar to previous example, but consider only values of 255 for filtering operation. Typical usage: dilate cloud values in input image that are flagged as 255
+
+\code
+pkfilter -i input.tif -o filter_stdev.tif -dz 3 -f median
+\endcode
+Filtering in spectral/temporal domain. Calculate the median value for each pixel, calculated on a moving window of width 3 (-dz 3) over all input bands. The output raster dataset will contain as many bands as the input raster dataset.
+
+\code
+pkfilter -i input.tif -o filter_stdev.tif -dz 1 -f stdev
+\endcode
+"Filtering" in spectral/temporal domain. No moving window (-dz 1). Calculate the standard deviation for each pixel, calculated on all input bands. The output raster dataset will contain a single band only.
+
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index 9894b39..b843991 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -62,7 +62,7 @@ int main(int argc, char *argv[])
Optionpk<string> label_opt("cn", "cname", "Name of the class label in the output vector dataset", "label");
Optionpk<short> geo_opt("g", "geo", "Use geo coordinates (set to 0 to use image coordinates)", 1);
Optionpk<short> down_opt("down", "down", "Down sampling factor (for raster sample datasets only). Can be used to create grid points", 1);
- Optionpk<short> buffer_opt("buf", "buffer", "Buffer for calculating statistics for point features ", 3);
+ Optionpk<short> buffer_opt("buf", "buffer", "Buffer for calculating statistics for point features ");
Optionpk<bool> disc_opt("circ", "circular", "Use a circular disc kernel buffer (for vector point sample datasets only, use in combination with buffer option)", false);
Optionpk<short> verbose_opt("v", "verbose", "Verbose mode if > 0", 0);
@@ -137,9 +137,7 @@ int main(int argc, char *argv[])
nvalid[it]=0;
ninvalid[it]=0;
}
- short theDim=(ruleMap[rule_opt[0]]==rule::point)? 1 : buffer_opt[0];
- if(verbose_opt[0]>1)
- std::cout << "boundary: " << buffer_opt[0] << std::endl;
+
ImgReaderGdal imgReader;
if(image_opt.empty()){
std::cerr << "No image dataset provided (use option -i). Use --help for help information";
@@ -913,6 +911,17 @@ int main(int argc, char *argv[])
assert(poGeometry!=NULL);
try{
if(wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ){
+
+ if(!buffer_opt.size()){
+ if(ruleMap[rule_opt[0]]==rule::point)
+ buffer_opt.push_back(1);//default
+ else
+ buffer_opt.push_back(3);//default
+ }
+
+ if(verbose_opt[0]>1)
+ std::cout << "boundary: " << buffer_opt[0] << std::endl;
+
OGRPolygon writePolygon;
OGRLinearRing writeRing;
OGRPoint writeCentroidPoint;
@@ -936,10 +945,34 @@ int main(int argc, char *argv[])
j_centre=static_cast<int>(j_centre);
i_centre=static_cast<int>(i_centre);
- double uli=i_centre-theDim/2-0.5;
- double ulj=j_centre-theDim/2-0.5;
- double lri=i_centre+theDim/2+0.5;
- double lrj=j_centre+theDim/2+0.5;
+ double uli=i_centre-buffer_opt[0]/2;
+ double ulj=j_centre-buffer_opt[0]/2;
+ double lri=i_centre+buffer_opt[0]/2;
+ double lrj=j_centre+buffer_opt[0]/2;
+ // double uli=i_centre-buffer_opt[0]/2-0.5;
+ // double ulj=j_centre-buffer_opt[0]/2-0.5;
+ // double lri=i_centre+buffer_opt[0]/2+0.5;
+ // double lrj=j_centre+buffer_opt[0]/2+0.5;
+
+ //nearest neighbour
+ ulj=static_cast<int>(ulj);
+ uli=static_cast<int>(uli);
+ lrj=static_cast<int>(lrj);
+ lri=static_cast<int>(lri);
+
+ if((polygon_opt[0]&&ruleMap[rule_opt[0]]==rule::point)||(ruleMap[rule_opt[0]]==rule::centroid)){
+ uli=i_centre;
+ ulj=j_centre;
+ lri=i_centre;
+ lrj=j_centre;
+ }
+
+ //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,uly;
@@ -948,7 +981,7 @@ int main(int argc, char *argv[])
if(polygon_opt[0]){
if(disc_opt[0]){
double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
- double radius=theDim/2.0*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
+ double radius=buffer_opt[0]/2.0*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
unsigned short nstep = 25;
for(int i=0;i<nstep;++i){
OGRPoint aPoint;
@@ -982,24 +1015,6 @@ int main(int argc, char *argv[])
}
}
- if(ruleMap[rule_opt[0]]==rule::centroid){
- uli=i_centre;
- ulj=j_centre;
- lri=i_centre;
- lrj=j_centre;
- }
- //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;
-
int nPointWindow=0;//similar to nPointPolygon for polygons
if(polygon_opt[0]){
if(writeTest)
@@ -1054,8 +1069,8 @@ int main(int argc, char *argv[])
thePoint.setX(theX);
thePoint.setY(theY);
- if(disc_opt[0]&&theDim>1){
- double radius=theDim/2.0*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
+ if(disc_opt[0]&&buffer_opt[0]>1){
+ double radius=buffer_opt[0]/2.0*sqrt(imgReader.getDeltaX()*imgReader.getDeltaY());
if((theX-x)*(theX-x)+(theY-y)*(theY-y)>radius*radius)
continue;
}
@@ -1198,12 +1213,15 @@ int main(int argc, char *argv[])
std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
if(class_opt.empty()){
- if(ruleMap[rule_opt[0]]==rule::point){//value at each point (or at centroid of polygon if line is set)
+ if(ruleMap[rule_opt[0]]==rule::point){//value at centroid of polygon
if(verbose_opt[0])
std::cout << "number of points in window: " << nPointWindow << std::endl;
for(int index=0;index<windowValues.size();++index){
//test
- assert(windowValues[index].size()==1);
+ if(windowValues[index].size()!=1){
+ cerr << "Error: windowValues[index].size()=" << windowValues[index].size() << endl;
+ assert(windowValues[index].size()==1);
+ }
double theValue=windowValues[index].back();
if(verbose_opt[0])
@@ -1235,6 +1253,9 @@ int main(int argc, char *argv[])
else{//ruleMap[rule_opt[0]] is not rule::point
double theValue=0;
for(int index=0;index<windowValues.size();++index){
+ if(windowValues[index].size()!=buffer_opt[0]*buffer_opt[0]){
+ cerr << "Error: windowValues[index].size()=" << windowValues[index].size() << endl;
+ }
if(ruleMap[rule_opt[0]]==rule::mean)
theValue=stat.mean(windowValues[index]);
else if(ruleMap[rule_opt[0]]==rule::stdev)
@@ -1659,7 +1680,7 @@ int main(int argc, char *argv[])
std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
if(class_opt.empty()){
- if(ruleMap[rule_opt[0]]==rule::point){//value at each point (or at centroid of polygon if line is set)
+ if(ruleMap[rule_opt[0]]==rule::point){//value at centroid of polygon
if(verbose_opt[0])
std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
for(int index=0;index<polyValues.size();++index){
@@ -2128,7 +2149,7 @@ int main(int argc, char *argv[])
std::cout << "write feature has " << writeCentroidFeature->GetFieldCount() << " fields" << std::endl;
}
if(class_opt.empty()){
- if(ruleMap[rule_opt[0]]==rule::point){//value at each point (or at centroid of polygon if line is set)
+ if(ruleMap[rule_opt[0]]==rule::point){//value at centroid of polygon
if(verbose_opt[0])
std::cout << "number of points in polygon: " << nPointPolygon << std::endl;
for(int index=0;index<polyValues.size();++index){
--
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