[R-sig-Geo] Problem in converting SpatialPolygonsDataFrame to owin object

[Roger Bivand]

On Fri, 15 Sep 2006, Adrian Baddeley wrote:

> 
> Debarchana Ghosh reports a problem in converting an object
> of class SpatialPolygonsDataFrame to class 'owin'.
> 
> Without any other information (without the data and knowing 
> nothing about 'sp') my guess is that you are trying to 
> make a single polygonal region out of several polygons that are adjoining 
> i.e. there are two polygons that share an edge, like the 
> borders of adjoining countries.

Yes, correct. The real culprit is the very messy shapefile (which was made 
available to me offline), which not only appears to be of 7 contiguous 
counties, but when dissolved using:

library(spgpc)
tcma_outer <- unionSpatialPolygons(as(tcma, "SpatialPolygons"), rep(1, 7))

has a single outer boundary and as many as 428 internal slivers. My 
immediate suggestion to the questioner would be to locate a better data 
source, such as:

http://www.census.gov/geo/cob/bdy/co/co00shp/co27_d00_shp.zip

In R:

library(rgdal)
tcma <- readOGR(".", "DG_tcma")
library(maptools)
mn <- readShapePoly("co27_d00.shp", proj4string=CRS("+proj=longlat +datum=WGS84"))
mn_FIPS <- paste(as.character(mn$STATE), as.character(mn$COUNTY), sep="")
match(as.character(tcma$COUNTY), mn_FIPS)
mn_subset <- mn[match(as.character(tcma$COUNTY), mn_FIPS),]
mn_subset_utm15N <- spTransform(mn_subset, CRS(proj4string(tcma)))
plot(mn_subset_utm15N, lwd=3, axes=TRUE)
plot(tcma, add=TRUE, border="red")
library(spgpc)
tcma_outer <- unionSpatialPolygons(as(mn_subset_utm15N, "SpatialPolygons"), rep(1, 7))
plot(tcma_outer)
library(spspatstat)
tcma_owin <- as(as(tcma_outer, "SpatialPolygons"), "owin")
plot(tcma_owin)

The key problem was the very messy shapefile, followed by not having 
understood that owin objects cannot be contiguous. The wrapper packages 
are on the sourceforge site. 

Hope this helps,

Roger

> 
> Did you want to end up with a single region of space,
> or with a list of different regions?
> 
> In the spatstat package, an object of class `owin' always represents 
> a region of space with an inside, an outside and a boundary. 
> To create a region with a polygonal boundary, you have to supply
> the coordinates of the *boundary* line segments. That is, every 
> edge of your polygon must be part of the *boundary* of the spatial region
> (the region separating the inside from the outside).
> 
> For example, to define the USA as an owin object, 
> you would supply the coordinates of the international borders
> (USA-Canada and USA-Mexico) and the ocean coast, but **not** the 
> interstate borders.
> 
> When you create a polygonal region as an object of class 'owin',
> spatstat checks its validity by checking for
> 	 - self-intersection in each polygon
> 	 - intersections between different polygons
> 	   (including common boundaries).
> It will not automatically `join' together two adjacent polygons.
> That kind of functionality is better handled by other code like gpclib.
> 
> Alternatively you can convert each separate polygon
> to an `owin' object, then use 'union.owin' to compute their union
> (which is doen by converting the polygons to binary images).
> 
> Adrian Baddeley
> 

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