[R-sig-Geo] Clipping grid using shapefile?

[Christopher Swingley]

Edzer,

* Edzer Pebesma <edzer.pebesma at uni-muenster.de> [2010-Mar-01 10:35 AKST]:
> >Reading the data in:
> >
> >  > library(rgdal)
> >  > library(sp)
> >  > vector <- readOGR("shapefiles", "watersheds")
> >  > ppt01 <- readGDAL("ppt01.asc")
> >  > proj4string(ppt01) <- CRS("+proj=aea +lat_1=55 +lat_2=65 +lat_0=50
> >    +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
> >    +towgs84=0,0,0")
>
> please tell us what was unexpected after you overlayed vector and
> ppt01.

When I loaded vector R said:

    with 23 features and 10 fields

When I overlay:

    > clip <- overlay(ppt01, vector)
    > class(clip)
    [1] "integer"
    > summary(clip)
    Min.  1st Qu.  Median   Mean   3rd Qu.  Max.
   1.000    5.000   7.000  7.564  11.000   23.000
    NA's
    20196990.000

I guess I was thinking I'd get something that had 23 rows in it, one for
each polygon in the vector file.

I tried something like this:

    > clip <- overlay(ppt01, vector[vector$SITECODE=='A-001',])

and got:

    > summary(clip)
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's
      1        1        1        1        1        1   20199300
    > mean(clip, na.rm = TRUE)
     [1] 1
    > clip[!is.na(clip)]
     [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


There's a total of (3505*5763 = 20199315 rows), so it's calculating it's
mean from 15 points.  When I use GDAL tools (gdal_translate to project
the ASC file to a TIF, gdal_rasterize to burn the shapefile polygon), I
also get 15 grid values, but the range of numbers is between 1691 to
2645 (not 1 to 1).

It seems like I'm very close, but am missing something simple.

Thanks,

Chris
-- 
Christopher S. Swingley
http://swingleydev.com/
<cswingle at gmail.com>