[R-sig-Geo] The length of common boundaries
Roger.Bivand at nhh.no
Thu Oct 13 18:39:15 CEST 2005
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  to generate this information (should be
> > the combination of the options 'sides' and 'length'.
> > Markus
> >  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
> 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
> 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 ...
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.
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
More information about the R-sig-Geo