[R-sig-Geo] The length of common boundaries
neteler at itc.it
Thu Oct 13 17:34:06 CEST 2005
(cc grass user list)
On Wed, Oct 12, 2005 at 05:43:35PM +0200, Markus Neteler wrote:
> On Tue, Oct 11, 2005 at 03:50:05PM +0200, Roger Bivand wrote:
> > On Tue, 11 Oct 2005, Yilin Liu wrote:
> > > Dear all,
> > >
> > > Any way to get the length of common boundaries shared by neighbouring
> > > counties in a map in R?
> > >
> > The short answer is no. The R polygon objects do not have topology.
> > The longer answer is that if you have e00 or ArcInfo binary vector layers,
> > then the polygons are built from lists of directed arcs, and the arcs have
> > lengths. So if your input data are in this format, it is feasible but is
> > not implemented.
> > It looks as though the GRASS6 vector format also provides similar
> > information. But topology is a GIS thing rather than an R thing, so no
> > solution within R is available.
> > If you use GRASS, it would be possible to look at how this might be done
> > (read a shapefile into GRASS, output an ASCII dump of the arcs, which
> > polygons they separate, and how long they are, and do something with this,
> > not losing trace of which polygon(s) belong to which observation - this
> > was where progress stopped the last time I looked). Arc lengths would be
> > nice for cartograms too.
> In GRASS 6 you can use v.to.db  to generate this information (should be
> the combination of the options 'sides' and 'length'.
>  http://grass.itc.it/grass60/manuals/html60_user/v.to.db.html
With the help of Radim I managed to find out how to do it.
It's pretty easy (but you need the today's version of GRASS 6.1-CVS
due to a related bugfix in v.db.addtable).
EXERCISE: HOW LONG ARE COMMON BOUNDARIES OF POLYGONS?
# Requires: GRASS 6.1-CVS from 13 Oct 2005 or later
# data: sudden infant deaths data from North Carolina
# data imported from SHAPE file with v.in.ogr
#let's have a look
#we work on a copy:
#we add a second layer to the map which references the boundaries of
#polygons. In the vector geometry we generate an ID (category) for each
v.category sids_nc out=sids_nc2 layer=2 type=boundary option=add
#we'll fetch the IDs (categories) of the polygons left and right from
#each boundary and store it into the attribute table linked to layer 2.
# cat_of_boundary | cat_of_left_polygon | cat_of_right_polygon | length_of_boundary
#We want only one category per boundary, that's why the sides check is
#needed (a boundary may consist of several pieces)
#So we create a new attribute table and link it to the new layer 2
#of the vector map:
v.db.addtable sids_nc2 layer=2 col="left integer,right integer,length integer"
#Now we query the polygon/boundary relationsships and store it into
#the attribute table linked to layer 2:
v.to.db map=sids_nc2 option=sides col=left,right layer=2
#Now we have unique categories for the boundaries and can calculate the
v.to.db map=sids_nc2 option=length col=length layer=2
#See the new attribute table containing the boundary lengths:
v.db.select sids_nc2 layer=2
# verification (let's check boundary #193):
d.vect sids_nc2 cat=193 layer=2 col=red type=boundary
# LEN: 12756.00 meters
#what does the attribute table say:
v.db.select sids_nc2 layer=2 | grep '^193'
#This is reasonably close since on screen digitization in d.measure
#isn't always that precise ...
More information about the R-sig-Geo