[R-sig-Geo] Map2poly and poly2nb with unconnected polygon areas

[Roger Bivand]

On Fri, 1 Dec 2006, Berta wrote:

> Hi R-sig-geo list,
> 
> A very easy question: I have a shapefile with 730 areas. Somo of these
> areas have two or three unconnected polygons, so that

Not so easy. There are 758 objects present, but some of the 730 
observations are associated with more than one object.

If you say (using the maptools package):

polys_758 <- readShapePoly("myfile.shp")

or using rgdal:

polys_758 <- readOGR("mydirectory", "myfile")

you will get a SpatialPolygonsDataFrame

Then (using maptools, and after installing gpclib), say:

polys_730 <- unionSpatialPolygons(as(polys_758, "SpatialPolygons"), 
  as.character(polys_758$cod))

you have 730 objects. The final question is whether the rows in the data 
frame are simply duplicated for the same cod value, and can be discarded, 
or whether they should be aggregated (sum, for example). 

df_758 <- as(polys_758, "data.frame")

convert to: df_730, set row.names(df_730) <- as.character(df_730$cod)

and:

out_730 <- SpatialPolygonsDataFrame(polys_730, data=df_730)

which will use the row names of df_730 to match against the 
dissolving/union criterion.

Then (using spdep):

x.nb <- poly2nb(as(out_730, "SpatialPolygons"))

should just work (and should put the cod values into the neighbour list 
IDs). Check with:

plot(out_730, border="grey")
plot(x.nb, coordinates(out_730), pch=".", col="blue")

It is quite typical for shapefiles to contain several objects for one 
observation, because the writing application is thinking in terms of 
geometry, not observations. Like having three of me, one for each 
telephone number, rather than one of me with three numbers, in a telephone 
number list.

Hope this helps,

Roger

> 
> length(myshapefile$att.data$cod) is 758 eventhough there are only 730 cods.
> 
> Now I want to obtain the neighbourhood structure using
> 
> x.2<-Map2poly(myshapefile, region.id="cod")
> x.nb<-poly2nb(x.2)
> sapply(x.nb,length) 
> 
> But it seems that it does not recognize the polygons that belong to the
> same area (indicated by cod), and again length(sapply(x.nb,length)) is
> 758 in stead of 730 . What am I doing wrong? Thanks a lot in advance,
> 
> Berta.
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Roger Bivand
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