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

Christopher Swingley cswingle at gmail.com
Tue Mar 2 17:24:47 CET 2010


All,

Thanks to a lot of help off-list, I managed to get it to work.  What I
was trying to do was clip an ASCII grid of climatological data to a
series of watershed polygons (stored in a shapefile) and calculate the
mean value for the watershed.

Here's the code that worked:

    > library(rgdal)
    > library(sp)
    > vector <- readOGR("Watersheds", "DENA")
    > ppt01 <- readGDAL("Climate/ppt01.asc")
    > clip <- overlay(ppt01, vector[vector$SITECODE=='DENA-001',])
    > mean(ppt01 at data[[1]][!is.na(clip)])
    > sd(ppt01 at data[[1]][!is.na(clip)])

I think my primary problem is not quite grasping how to access the
various parts of the data objects loaded.  More reading of Data
Manipulation with R and Applied Spatial Data Analysis with R.

I may also try an approach using the raster library to see which is
faster for my data.  A python/GDAL version I wrote takes about 6 seconds
for each grid x polygon intersection calculation.

Cheers,

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



More information about the R-sig-Geo mailing list