[R-sig-Geo] The length of common boundaries

[Markus Neteler]

(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 [1] to generate this information (should be
> the combination of the options 'sides' and 'length'.
> 
> Markus
> 
> [1] 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
 d.mon x0
 d.vect sids

#we work on a copy:
 g.copy vect=sids,sids_nc

#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
#boundary:
 v.category sids_nc out=sids_nc2 layer=2 type=boundary option=add

#Underlying idea:
#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.
#In general:
# 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
#lengths:
 v.to.db map=sids_nc2 option=length col=length layer=2 

#Done.

#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
 d.zoom
 d.measure
# LEN:     12756.00 meters

#what does the attribute table say:
 v.db.select sids_nc2 layer=2 | grep '^193'
#190|65|68|12814

#This is reasonably close since on screen digitization in d.measure
#isn't always that precise ...


Best

 Markus