[R-sig-Geo] Calculating lengths of shared edges for neighbouring polygons (in R, without resorting to GRASS)

[Graeme Cumming]

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