[GRASS-stats] Re: Problem with v.distance in spgrass6 package in R
Roger Bivand
Roger.Bivand at nhh.no
Thu Nov 3 03:11:09 EDT 2011
On Wed, 2 Nov 2011, toke wrote:
> Hi Roger and everybody else
>
> I finally found time to translate the problem into spearfish data.
>
> To recapitulate I would like to use the v.distance command to transfer
> information from polygons to point data.
>
> The problem is that the v.distance command do not update the column
> specified to receive the calculation done by the v.distance command.
>
I do not see the same problem, and the command works for me - can others
please try to reproduce it?
Roger
My output in R:
> library(spgrass6)
Loading required package: sp
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.8.1, released 2011/07/09
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files: (autodetected)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: GRASS 6.4.2svn47586 (2011)
and location: spearfish60
> execGRASS("g.copy", vect="'fields at PERMANENT',fields1")
Copy vector <fields at PERMANENT> to current mapset as <fields1>
> execGRASS("g.copy", vect="'archsites at PERMANENT',archsites1")
Copy vector <archsites at PERMANENT> to current mapset as <archsites1>
> execGRASS("v.db.addcol", map="archsites1",
+ columns="\"fields double precision\"")
> fields1<-readVECT6("fields1")
WARNING: The map contains islands. To preserve them in the output map,
use
the -c flag
Exporting 65 areas (may take some time)...
100%
v.out.ogr complete. 67 features written to <fields1> (ESRI_Shapefile).
OGR data source with driver: ESRI Shapefile
Source: "/home/rsb/topics/grassdata/spearfish60/rsb/.tmp/reclus2", layer:
"fields1"
with 67 features and 2 fields
Feature type: wkbPolygon with 2 dimensions
> archsites1<-readVECT6("archsites1")
Exporting 25 geometries...
100%
v.out.ogr complete. 25 features written to <archsite> (ESRI_Shapefile).
OGR data source with driver: ESRI Shapefile
Source: "/home/rsb/topics/grassdata/spearfish60/rsb/.tmp/reclus2", layer:
"archsite"
with 25 features and 3 fields
Feature type: wkbPoint with 2 dimensions
> summary(archsites1)
Object of class SpatialPointsDataFrame
Coordinates:
min max
coords.x1 589860 608355
coords.x2 4914479 4926490
Is projected: TRUE
proj4string :
[+proj=utm +zone=13 +datum=NAD27 +units=m +no_defs +ellps=clrk66
+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat]
Number of points: 25
Data attributes:
cat str1 fields
Min. : 1 No Name :12 Min. : NA
1st Qu.: 7 Bob Miller : 1 1st Qu.: NA
Median :13 Boulder Creek Cabin: 1 Median : NA
Mean :13 Canyon Station : 1 Mean :NaN
3rd Qu.:19 Cole Creek Mine : 1 3rd Qu.: NA
Max. :25 Elkhorn Peak : 1 Max. : NA
(Other) : 8 NA's : 25
> set.echoCmdOption(TRUE)
[1] FALSE
> execGRASS("v.distance", flags="p", from="archsites1", to="fields1",
from_type="point", to_type="point,line,area", dmax=as.integer(1),
upload="to_attr", column="fields", to_column="cat")
GRASS command: v.distance -p from=archsites1 to=fields1 from_type=point
to_type=point,line,area dmax=1 upload=to_attr column=fields to_column=cat
100%
100%
100%
from_cat|fields
1|63
2|63
3|63
4|63
5|null
6|null
7|null
8|null
9|25
10|63
11|null
12|null
13|63
14|63
15|null
16|63
17|null
18|63
19|63
20|63
21|null
22|63
23|63
24|null
25|63
v.distance complete.
> set.echoCmdOption(FALSE)
[1] TRUE
> archsites1<-readVECT6("archsites1")
Exporting 25 geometries...
100%
v.out.ogr complete. 25 features written to <archsite> (ESRI_Shapefile).
OGR data source with driver: ESRI Shapefile
Source: "/home/rsb/topics/grassdata/spearfish60/rsb/.tmp/reclus2", layer:
"archsite"
with 25 features and 3 fields
Feature type: wkbPoint with 2 dimensions
> summary(archsites1)
Object of class SpatialPointsDataFrame
Coordinates:
min max
coords.x1 589860 608355
coords.x2 4914479 4926490
Is projected: TRUE
proj4string :
[+proj=utm +zone=13 +datum=NAD27 +units=m +no_defs +ellps=clrk66
+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat]
Number of points: 25
Data attributes:
cat str1 fields
Min. : 1 No Name :12 Min. :25.00
1st Qu.: 7 Bob Miller : 1 1st Qu.:63.00
Median :13 Boulder Creek Cabin: 1 Median :63.00
Mean :13 Canyon Station : 1 Mean :60.47
3rd Qu.:19 Cole Creek Mine : 1 3rd Qu.:63.00
Max. :25 Elkhorn Peak : 1 Max. :63.00
(Other) : 8 NA's :10.00
> table(archsites1$fields)
25 63
1 14
## alternatively:
> execGRASS("v.distance", from="archsites1", to="fields1",
from_type="point", to_type="point,line,area", dmax=as.integer(1),
upload="to_attr", column="fields", to_column="cat")
GRASS command: v.distance from=archsites1 to=fields1 from_type=point
to_type=point,line,area dmax=1 upload=to_attr column=fields to_column=cat
100%
100%
100%
25 categories read from the map
25 categories exist in the table
25 categories read from the map exist in the table
25 records updated
v.distance complete.
## also OK (overwriting field values)
> sessionInfo()
R version 2.14.0 (2011-10-31)
Platform: x86_64-unknown-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spgrass6_0.7-4 XML_3.4-3 rgdal_0.7-1 sp_0.9-91
loaded via a namespace (and not attached):
[1] grid_2.14.0 lattice_0.20-0
> Here is my script
>
> grass
>
> g.gisenv
> g.mapset mapset=user1
>
> g.copy vect='fields at PERMANENT',fields1 # adding polygon to mapset
>
> g.copy vect='archsites at PERMANENT',archsites1 # adding points to mapset
>
> v.db.addcol map=archsites1 columns="fields double precision" # adding
> column to points data
>
> v.distance from=archsites1 to=fields1 from_type=point
> to_type=point,line,area dmax=1 upload=to_attr column=fields to_column=cat #
> adding cat information from polygon to point data
>
> ### v.distance works in GRASS ###
>
> ### trying to do the same thing in R ###
>
> v.db.dropcol archsites1 column=fields ## removing column
>
> v.db.addcol map=archsites1 columns="fields double precision" ## adding a
> clean column
>
>
> R # starting R within grass
>
> library(spgrass6) # calling library
>
> gmeta6()
>
> fields1<-readVECT6("fields1") # reading polygon into R
>
> archsites1<-readVECT6("archsites1") # reading points into R
>
> archsites1 # looking at data
>
> # Trying to do the same v.distance calculation in R.
> ### The specified column is not updated ###
>
> execGRASS("v.distance", from="archsites1",
> to="fields1", from_type="point", to_type="point,line,area",
> dmax=as.integer(1),
> upload="to_attr", column="fields", to_column="cat")
>
>
> ### However if the flags “p” is add to the v.distance calculation it is
> possible to see that v.distance
> ### actually do the calculations.
>
> execGRASS("v.distance", flags="p", from="archsites1",
> to="fields1", from_type="point", to_type="point,line,area",
> dmax=as.integer(1),
> upload="to_attr", column="fields", to_column="cat")
>
> ### script finished ###
>
> There seems to be a bug in spgrass6.
>
> Cheers Toke
>
>
> --
> View this message in context: http://osgeo-org.1803224.n2.nabble.com/Problem-with-v-distance-in-spgrass6-package-in-R-tp6939286p6955833.html
> Sent from the Grass - Stats mailing list archive at Nabble.com.
> _______________________________________________
> grass-stats mailing list
> grass-stats at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-stats
>
--
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the grass-stats
mailing list