[R-sig-Geo] extracting raster data with point and polygons

[Oscar Perpiñan]

Hello,

If I understood you right, I think you can solve it using `zonal` and
`mask`. This is my proposal:

library(raster)

###### Your code
#define DEM
dem2m=raster(x=matrix(runif(22500,min=2000,max=2500),nrow=150),crs='+proj=utm
+zone=13 +datum=NAD83')
extent(dem2m)=c(0,300,0,3000)
#
#define satellite image
fsca.dem=raster(x=matrix(runif(100,0,100),nrow=10),crs='+proj=utm
+zone=13 +datum=NAD83')
extent(fsca.dem)=c(0,300,0,3000)

#define point data
obs=data.frame(x=runif(5,0,300),y=runif(5,0,300),sd=runif(5,0,5))
coordinates(obs)=~x+y
proj4string(obs)='+proj=utm +zone=13 +datum=NAD83'

###### My proposal
## Create a Raster with the same resolution as `dem2m`, whose values are
## the cell numbers of `fsca.dem`
xyDEM2m <- xyFromCell(dem2m, cell = seq_len(ncell(dem2m)))
zones <- setValues(dem2m, cellFromXY(fsca.dem, xyDEM2m))
## Compute the cv for each zone. `zoneVariation` has the same
## resolution as `fsca.dem`
zoneVariation <- setValues(fsca.dem, zonal(dem2m, zones, fun = cv)[,2])
## We are only interested in those cells with observations
varObs <- mask(zoneVariation, obs)

If you work with large rasters maybe it is better to apply mask before.

Best,

Oscar.
-----------------------------------------------------------------
Oscar Perpiñán Lamigueiro
Dpto. Ingeniería Eléctrica (ETSIDI-UPM)
Grupo de Sistemas Fotovoltaicos (IES-UPM)
URL: http://oscarperpinan.github.io
Twitter: @oscarperpinan


2014-07-22 2:26 GMT+02:00 Dominik Schneider <Dominik.Schneider at colorado.edu>:
> Hi -
> I have a 2 meter DEM, a 30m satellite image and some point data. Basically,
> I would like to extract the 2m dem pixels that correspond with the 30m
> satellite data, but only for the satellite pixels that correspond with a
> point observation. i've been trying to get at the answer from the sp
> vignette 'over' but am stuck.
> an example:
> #define DEM
> dem2m=raster(x=matrix(runif(22500,min=2000,max=2500),nrow=150),crs='+proj=utm
> +zone=13 +datum=NAD83')
> extent(dem2m)=c(0,300,0,3000)
> #
> #define satellite image
> fsca.dem=raster(x=matrix(runiff(100,0,100),nrow=10),crs='+proj=utm +zone=13
> +datum=NAD83')
> extent(dem2m)=c(0,300,0,3000)
> #
> #define point data
> obs=data.frame(x=runif(5,0,300),y=runif(5,0,300),sd=runif(5,0,5))
> coordinates(obs)=~x+y
> proj4string(obs)='+proj=utm +zone=13 +datum=NAD83'
> #
> #--- in my case, I resampled the sat image to the dem grid. but thats
> already taken care of in the example.
> ##dem.pe=projectExtent(dem2m,crs=projection(dem2m))
> ##res(dem.pe) <- 30
> ##fsca.dem=projectRaster(fsca,dem.pe,method='ngb')
> #
> # This is where I am lost
> # ----------- I thought I'd use polygons to extract the dem data. The catch
> is, rather than extracting every polygon of data (which takes longer than
> I've been willing to wait given 2e6 pixels) I would like to subset this
> spatialpolygondf to the 5 polygons that have observational data. if the
> same polygon has multiple obs in it, it should be replicated.
> fsca.poly=rasterToPolygons(fsca.dem)
>
> #first i tried this but this extracts the fsca values. I need this too
> eventually but was hoping for a subsetted spatialpolygon.dataframe
> fsca.pts=over(ss.depth,fsca.poly)
>
> #second i tried this but the list has rownames that I can't interpret.
> (fsca cellnumbers?)
> fsca.pts.poly=over(ss.depth,fsca.poly,returnList=T)
> polynames=unlist(lapply(fsca.pts.poly,rownames))
>
> #but what i'd like to do is apply the cv function on the 225 2m elev pixels
> that correspond to the 30 fsca pixel for which i have an observation point.
> elev.cv=extract(dem2m,*y=fsca.poly.w.pts*,fun=cv)
>
> # as a note, the buffer option in extract would not work currently because
> the point data isn't necessarily in the middle dem pixel relative to the
> fsca pixel. so one solution might be to 1. subset the polygon dataframe 2.
> use the cell center from the fsca pixel to extract from the dem raster with
> a buffer of 7 cells.
>
> anyway, i appreciate any help you can give, in particulary with subsetting
> the spatialpolygondf.
> Thanks
> Dominik
>
>
>> sessionInfo()
> R version 3.0.2 (2013-09-25)
> Platform: x86_64-apple-darwin10.8.0 (64-bit)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] maptools_0.8-30    reshape2_1.2.2     plyr_1.8
> rgdal_0.8-16
> [5] ggplot2_0.9.3.1    RColorBrewer_1.0-5 scales_0.2.3
> raster_2.2-31
> [9] sp_1.0-15
>
> loaded via a namespace (and not attached):
>  [1] MASS_7.3-29      colorspace_1.2-4 compiler_3.0.2   dichromat_2.0-0
>  [5] digest_0.6.4     foreign_0.8-55   grid_3.0.2       gtable_0.1.2
>  [9] labeling_0.2     lattice_0.20-29  munsell_0.4.2    proto_0.3-10
> [13] stringr_0.6.2    tcltk_3.0.2      tools_3.0.2
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo