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

[Roger Bivand]

On Tue, 30 Oct 2007, Agustin Lobo wrote:

> That's really helpful, thanks.
>
> For "everyday" use, I've made:
>
> "slotspol" <- function(pol, sl="ID")
> sapply(slot(pol, "polygons"), function(x) slot(x, sl))
>
> Which can be used for any slot name:
> slotspol(Up)
> slotspol(Up,"area")
>
> Also, few questions:
> 1. I observe that pols$Site is identical to pols at data$Site. Is not 
> pols at data$Site the correct syntax? I mean that we retrieve variable "Site" 
> from slot "data" in object "pols".

Both are possible, but "$" methods are provided for all Spatial*DataFrame 
objects so that they behave like regular data frames in many settings, so 
that coercing to data frame becomes unnecessary. Using the "@" directly is 
not necessarily safe, so in particular please use the "$<-" method for 
writing to (possibly new) variables in the Spatial*DataFrame. To see how 
"$" works, try using a Spatial*DataFrame as the data= argument in lm() or 
glm(), etc. It just "works".

>
> 2. In your code, I don't understand
> sites_Lupuna <- length(grep("Lupuna", as.character(pols$Site))) > 0
>
> because the result is just T (the result of length(grep("Lupuna", 
> as.character(pols$Site)))  is 3)
> and thus subsetting with sites_Lupuna yields an object that is identical
> to the original:

I don't have your data, so I cannot set up the condition vector exactly. 
Here I was trying to be too clever, you want a vector of TRUE or FALSE 
values of length 140, and this condition wasn't the right one.

cond <- "Aleman" %in% as.character(pols$Site)

should just have TRUE for Aleman, not otherwise. Your understanding below 
is OK.

Roger

>
>>  pols_Lupuna <- pols[sites_Lupuna,]
>>  length(pols_Lupuna at polygons)
> [1] 140
>>  length(pols at polygons)
> [1] 140
>
> Instead:
>>  sites_Lupuna <- grep("Lupuna", as.character(pols$Site))
>>  sites_Lupuna
> [1] 1 2 3
>>  pols_Lupuna <- pols[sites_Lupuna,]
>
> seems to provide a correctly subset of the original pols object ("Lupuna 1", 
> "Lupuna 2", "Lupuna 3")
>
> Am I wrong? I feel quite insecure on this.
>
> Agus
>
> Roger Bivand escribió:
>>  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