[R-sig-Geo] Subsetting SpatialPolygonsDataFrame and empty ID slots

[Roger Bivand]

On Mon, 29 Oct 2007, Agustin Lobo wrote:

> Dear list,
>
> I have imported a shp file with:
>
>> pols <- readOGR("../AllTransectPolygons02",
> layer="AllTransectPolygons02")
>
>> class(pols)
> [1] "SpatialPolygonsDataFrame"
> attr(,"package")
> [1] "sp"
>
>> pols at data
>     AREA PERIMETER               Site
> 0   2500      1010           Lupuna 1
> 1   2500      1010           Lupuna 3
> 2   2500      1010           Lupuna 2
> 3   2500      1010             Aleman
> 4   2500      1010        Zancudoyacu
> etc
>
> and now I want to select few polygons according
> to a vector of Site names.
>
> Thinking on how doing this, I thought I could select
> the corresponding IDs in the data table and then use
> the ID slot in the polygons, but
> I've found that the polygons themselves actually
> lack the ID information:
>

Please treat all the long named access thingies in sp as deprecated - they 
will be flagged as such in the next release, and removed one release after 
that. S4 is not like that.

library(sp)
library(rgdal)
example(readOGR)
class(Up)
getSlots(class(Up))
class(slot(Up, "polygons")[[1]])
getSlots(class(slot(Up, "polygons")[1]))

So:

sapply(slot(Up, "polygons"), function(x) slot(x, "ID"))

>> getPolygonsIDSlot(pols)
> Error in getPolygonsIDSlot(pols) :
>   no slot of name "ID" for this object of class "SpatialPolygonsDataFrame"
>
> How does sp relate the data table to the polygons if there is no ID?
>

Of course there is an ID, but you can get at either by:

sapply(slot(pols, "polygons"), function(x) slot(x, "ID"))

or

rownames(as(pols, "data.frame"))


> In any case, which would the most efficient way of subseting
> the SpatialPolygonsDataFrame into another SpatialPolygonsDataFrame
> with the selected polygons?

Returning to the original question, subsetting is simply by the "[" 
operator, just like any other data.frame:

sites_Lupuna <- length(grep("Lupuna", as.character(pols$Site))) > 0
summary(sites_Lupuna)
pols_Lupuna <- pols[sites_Lupuna,]

Just think of Spatial*DataFrame objects as data.frame objects and things 
will be much clearer.

Roger

> Should I convert to a maptools object
> and use subset.polylist()? If so, how do I convert? I've tried:
>
> pols.maptools <- SpatialPolygons2PolySet(pols)
>

PS. PolySet objects are really for drawing coastlines in this context.

> but the data table is lost.
>
> Thanks!
>
>

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