An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140721/07788041/attachment.pl>
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 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
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140722/9228da5f/attachment.pl>
Hi,
Could you clarify where you would use mask() to make this more efficient for larger rasters? It seems this will always use all the cellnumbers. I though maybe mask(zones,obs) but then zonal() threw an error.
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.