The length of common boundaries

(cc grass user list)

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