[R-sig-Geo] Calculating lengths of shared edges for neighbouring polygons (in R, without resorting to GRASS)
Roger Bivand
Roger.Bivand at nhh.no
Wed Jun 29 15:01:50 CEST 2011
On Wed, 29 Jun 2011, Graeme Cumming wrote:
> 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:
Could I suggest using ArcGIS - convert the vector feature class to a
vector coverage, then use the AVCBin driver in rgdal to read the PAL and
ARC versions? Refer to:
http://www.asdar-book.org/exercises.php?excode=3
from ASDAR, download http://www.asdar-book.org/bundles/nb_ARC.zip, unzip,
then srarting R in Arc/:
library(rgdal)
sy <- readOGR(dsn="syr_cov1", layer="PAL", drop_unsupported_fields=TRUE)
sy_arc <- readOGR(dsn="syr_cov1", layer="ARC",
drop_unsupported_fields=TRUE)
sy_arc$length <- SpatialLinesLengths(sy_arc)
summary(sy_arc)
gives the boundary lengths by right and left neighbour. You may need to
look at and modify coverage2nb.R to see how to get the correct FID values
- the function as it is only outputs neighbours. Note that *POLY_ 1 is the
excluded background. This may be simpler than trying to get the lengths in
R. The alternative is to use GRASS, as in the previous additional material
entry on the book website.
Hope this helps,
Roger
>
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
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
mailing list