[GRASSLIST:8605] Re: [R-sig-Geo] The length of common boundaries

[Roger Bivand]

On Thu, 13 Oct 2005, Markus Neteler wrote:

> (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 ...
> 

Apart from changing the declaration of length from integer to double, this 
really does the job very well. Thank you, Markus and Radim! The output of 
v.db.select sids_nc2 layer=2 gives all that is needed to construct the 
needed files, a *.gal contiguous neighbour file, a *.gwt file with 
neighbours and lengths per neighbour, and a simple perimeter file per 
unique ID. Something along these lines will reach spgrass6 soon.

Best wishes,

Roger


> 
> Best
> 
>  Markus
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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