Skip to content

extracting raster data with point and polygons

4 messages · Dominik Schneider, 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 mean that you may save computation time if you apply `cv` (or
whatever function) only over those cells that correspond to the point
observations. However, this does not work with `zonal` as you point.

Best,

Oscar.