[R-sig-Geo] Problem with poly2nb after unionSpatialPoly

James Rooney ROONEYJ4 at tcd.ie
Mon Jul 8 18:11:04 CEST 2013


Thanks Roger - I will take a look at chapter 5 and spCbind(). I figured I was scarmbling the order somewhere along the way but could not figure where exactly.
Many thanks for your help.

James

________________________________________
From: Roger Bivand [Roger.Bivand at nhh.no]
Sent: 08 July 2013 16:46
To: James Rooney
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Problem with poly2nb after unionSpatialPoly

On Mon, 8 Jul 2013, Roger Bivand wrote:

> On Mon, 8 Jul 2013, James Rooney wrote:
>
>> Hi all,
>>
>> I am trying to setup a geospatial CAR model to implement areal spatial
>> smoothing (similar to Fig11.11 in "Applied Spatial Data Analysis with R"
>> except without the beta term).
>> I am using ESRI shapfiles downloaded from Census Ireland
>> (http://www.cso.ie/en/census/census2011boundaryfiles/)
>> As some of the boundaries have changed over various census years, I need to
>> join certain polygons - note these are adjacent polygons and I have used
>> unionSpatial Polygons to try to do this. I am using a join_list of polygons
>> that have changed and are to be joined - this is just a CSV file I made
>> myself. The problem I am having is that my code is not joining the correct
>> polygons and is in fact joining non-adjacent unrelated polygons. I can't
>> figure out how to fix this.
>>
>> My code looks like this: (FYI - EDdata contains a unique code for each
>> polygon in a factor variable called GEOGID).
>>
>> EDdata<-readOGR(".","Census2011_Electoral_Divisions_generalised20m",TRUE)
>> names(EDdata)
>> join_list<- read.csv("joins.csv", header=TRUE)
>>
>> ED2<-merge(EDdata,join_list,by="GEOGID",all=TRUE)

OK. Don't do this, use the maptools spCbind() method after first having
set the row.names of EDdata and join_list (which is a data.frame object).
They are almost certainly not in the right order. Check the order with
match(), and if necessary change the order of join_list to match exactly.

>
> What is merge() - which package? If you mean the merge method in sp, it is
> not intended for joining geometries, nor I think is the earlier method in
> raster.
>
> The standard approach since forever is unionSpatialPolygons() in maptools, or
> eqivalently gUnaryUnion() in rgeos. Unless this merge() function is making
> strong assumptions, it will mess up the aggregation of multiple observed
> values from the joined/unioned Polygons objects. There have been multiple
> threads on this list on the use of gUnaryUnion().
>
> Hope this clarifies,
>
> Roger
>
>>
>> # modify copied geogids as per join list
>> ED2$GEO<-as.character(ED2$GEOGID)
>> ED2$GEO[!is.na(ED2$Joined)]<-as.character(ED2$Joined[!is.na(ED2$Joined)])
>> ED2$GEOGID2<-as.factor(ED2$GEO)
>>
>> # apply join function to polygons
>> EDdata$GEOGID2<-ED2$GEOGID2
>> EDpoly_joined <- unionSpatialPolygons(EDdata, EDdata$GEOGID2)
>>

At this point EDdata$GEOGID2 is a factor, and is in the wrong order. If
you use row.names() to set and retreive the character versions, it will be
easier. Maybe the new vignette in maptools (two sections of ch. 5 in the
book you refer to) would help, or since you have the book, you could read
the sections there.

Roger

>>
>> When I plot specific joined polygons they are not the adjacent ones I
>> intended to join. As far as I can see - GEOGID and GEOGID2 match each other
>> so I am a loss as to why the unionSpatialPolygons function appears to join
>> to incorrect polygons.
>>
>> FYI the join list looks like this (its longer than this but this should
>> give you an idea).
>>
>>                               Joined            variable       GEOGID
>> E01001                   E01001/01019     Sub1       E01001
>> E10001                         E10001         Sub1       E10001
>> E10004                         E10004         Sub1       E10004
>> E11001                   E11001/11026     Sub1       E11001
>> E11002                   E11002/11055     Sub1       E11002
>> E11003                   E11003/11092     Sub1       E11003
>>
>>
>>
>> I cannot see where I am going wrong -I'd be very grateful for any
>> assistance/advice.
>>
>> Many thanks,
>>
>> James
>> _______________________________________________
>> 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