[R-sig-Geo] Calculating lengths of shared edges for neighbouring polygons (in R, without resorting to GRASS)
Graeme Cumming
gscumming at googlemail.com
Wed Jun 29 13:41:28 CEST 2011
Dear all,
I am trying to calculate lengths of shared edge for neighbouring
polygons in a shapefile. I have most of the code working, I think, but
I don't get the right answer (estimated edge lengths are too long by
about half as much again). I'd be grateful if someone can help me
figure out where I'm going wrong.
To offer some background: My underlying motive is to generate a
boundary length file for use with the UQ MARXAN conservation planning
software. I have found a number of related queries on this list and as
best I can tell, the only viable solution to this problem so far
proposed has been to go via GRASS. I am teaching my students
conservation planning in ArcGIS, however, and I don't want them to
have to grapple with other GIS packages just yet (and yes, I am aware
of the various other alternatives for boundary file creation - most of
them don't work in ArcGIS 9.3 or 10, ArcView won't run in Windows 7,
and Zonae Cogito represents too niche-oriented a step away from
classical GIS for my class).
Each row in the output file needs to contain IDs of polygons 1 and 2,
which are neighbours, followed by the length of shared boundary. The
input shapefie is currently a vector grid of squares generated in
ArcGIS using fishnet or Hawth's tools. The grid is clipped by the
boundary of the western cape, so there is a set of irregular edge
cells, but central cells are square. I am testing this with a data set
in Albers Equal area projection and a cell side length of 10,000m. I
should also add that I am working on a Mac, so rgeos is not currently
available.
The neighbour list element of the problem seems to be working fine. I
have coded it using poly2nb as follows:
library(spdep)
library(gpclib)
library(sp)
library(PBSmapping)
#call shapefile and plot it to check that it is functional
wcape <-readShapePoly("WCapeGrid1_Clip"); #then define the projection
proj4string(wcape)<-CRS("+proj=aea +datum=WGS84 +units=m");
plot(wcape);#check that the map looks OK
nbs<-poly2nb(wcape, queen=FALSE); #we want shared edges, not just
vertices
#use 'summary(nbs)' to summarize contents of nbs
wcapenlist<-nb2listw(nbs,style='B');
#use 'names(wcapenlist)' to view variables in the data structure
nbs2<-wcapenlist$neighbours #extracts neighbour data from object
#then we do some fiddling around to get the data in more workable
format:
nbs3 <- nb2listw(nbs2, style="B") #generate a list
nbs4 <- as(as_dgRMatrix_listw(nbs3), "CsparseMatrix") #coerce it into
an nrowsxnrows sparse matrix
nbs5<-which(nbs4!=0,arr.ind=T)#find and make a list of row, col IDs of
nonzero entries
nrows <- nrow(nbs5); #find how long the matrix is
duplicates <- matrix(data=0,nrow=nrows); #create and add a colum to
mark duplicates
nbs6<-cbind(nbs5,duplicates);
#loop through the matrix to find complementary pairs (a,b and b,a) and
mark them with '1'
for (j in 1:nrows) {
for (k in j:nrows) {
if (j!=k && nbs6[j,1]==nbs5[k,2] && nbs6[j,2]==nbs6[k,1])
nbs6[k,3]<- 1;
}
}
#create a new matrix that contains only unique pairs
#this step can probably be done more efficiently
nuniquerows <- nrows-sum(nbs6[,3]);
nbs7 <- matrix(data=0,nrow=nuniquerows,ncol=2);
ticker<-0;
for (j in 1:nrows){
if(nbs6[j,3]==0){
ticker<-ticker+1;
nbs7[ticker,1]=nbs6[j,1];
nbs7[ticker,2]=nbs6[j,2];
}
}
OK. The problem comes somewhere in the next set of steps:
#we now have a list of all neighbour IDs. The next step is to
calculate shared boundary lengths.
#we do this by (1) picking pairs to compare; (2) unioning pairs of
polygons, which effectively removes the shared boundary; and then (3)
subtracting the perimeter of the union from the summed perimeters of
the two individual polygons.
###############################################################################################
nbscommon<- matrix(data=0,nrow=nuniquerows,ncol=1) #create a vector to
hold the lengths of common (shared) edges
for (j in 1:nuniquerows){
poly1=wcape[nbs7[j,1],];#extract first polygon
poly2=wcape[nbs7[j,2],];#extract its neighbour
poly1PBS <- combinePolys(SpatialPolygons2PolySet(poly1)); #convert to
PBS mapping format
poly2PBS <- combinePolys(SpatialPolygons2PolySet(poly2));
union.p1p2 <- combinePolys(joinPolys(poly1PBS, poly2PBS,"UNION",
maxVert = 1000));
#plotPolys(union.p1p2); addPolys(poly1PBS); #check that the union
looks OK
lengthp1 <- calcLength(poly1PBS); #calculate perimeter of poly1
lengthp1 <- sum(lengthp1$length); #use 'sum' because edges sometimes
lead to split polys
lengthp2 <- calcLength(poly2PBS); #calculate perimeter of poly2
lengthp2 <- sum(lengthp2$length);
lengthunion <- calcLength(union.p1p2); #calculate perimeter of union
of poly1 and poly2
lunion <- sum(lengthunion$length);
nbscommon[j] <- (lengthp1+lengthp2-lunion)/2; #shared/common perim is
half the difference between union and sum of individual polys
}
boundarymat <- cbind(nbs7,nbscommon);
save (boundarymat,file="boundary.txt",ascii=TRUE);
#export boundarymat to text file and we're done.
So far so good, I would have thought. But the problem is that my
square cells have side length 10km, and these cells are assigned a
side length of 15km by my algorithm. What aren't I seeing? Does
'combinePolys' double up on lines? I have a similar algorithm working
fine in Matlab, main problem being that my students don't have access
to the necessary toolboxes and so can't run it for themselves.
Thanks for your help!
Graeme Cumming
Prof. Graeme S. Cumming
Percy FitzPatrick Institute
University of Cape Town
Tel. +27-21-650-3439
http://www.fitzpatrick.uct.ac.za/gcumming/index.htm
More information about the R-sig-Geo
mailing list