What is the most efficient way to do a spatial join of all raster points within a polygon layer area?

Question

I have a raster GeoTIFF of a digital elevation model and a polygon shapefile with spatial boundaries of the raster. Here is a snapshot of raster and polygon layer:

The shapefile has one attribute named "plot #". My ultimate goal is to extract individual pixel values in the raster that fall under the polygons and join them onto the new shapefile that has plot# field. Currently, my workflow does this: In Qgis, export raster to point file (huge file) >> perform spatial join over the plot boundaries using custom programs in R. The final output table should look like (bonus if it is a shapefile!):

Plot# height

1 60

1 30

1 40

1 .

1 .

2 20

2 30

. .

Now the above workflow has a few issues. First, to me, the raster-to-point file export seems unnecessary and a slow process. Second, the spatial join is painfully slow for 1800 polygons on a dataframe. I wonder if: 1. there is any way to bypass the raster-to-point conversion step and directly apply the spatial join on raster pixels similar to what 'zonal statistics' plugin in this does but for all individual raster points within the polygon. I find zonal statistics to be very efficient and fast but it outputs only summary stats, not individual points. 2. If there is no other alternative, is it possible to do the same steps i.e. raster-to-point >> spatial join all in R instead of going back and forth in this and R?

How should I tackle the issue?

Answer 1

This can easily be accomplished in R using raster::extract or, for a much faster option, the velox package. There is no need to convert to an ASCII or point format, just read in your raster using the raster package (raster, stack or brick functions).

Here we create some example data and plot it.

    library(velox)
    library(raster)

    r <- raster(ncol=100, nrow=100)
      r[] <- runif(ncell(r))
    polys <- spPolygons(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)), 
                        rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)))
    polys <- SpatialPolygonsDataFrame(polys, data.frame(ID=1:2))    
    plot(r)
      plot(polys, add=TRUE)

Now, using raster::extract we can extract the raster values into a list object where the list elements represents the values in each polygon. Once the values are extracted we use lapply (list apply) to apply a function and return a statistic or results of a custom function. In this example we return the mean and then write a simple function that returns the number of values above a threshold.

( v <- extract(r, polys) )
unlist(lapply(v, FUN=mean))

f <- function(x, p=0.5) { length(x[x >= p]) } 
unlist(lapply(v, FUN=f))

We can join the results to the @data slot, holding the attribute data such as plot ID, of our polygons.

pmean <- unlist(lapply(v, FUN=mean)
polys@data <- data.frame(polys@data, mean = pmean))
polys@
spplot(polys, "mean")

Here is a quick velox example where we coerce our raster into a velox object and then return the mean.

vx <- velox(r)
( ex.mat <- vx$extract(polys, fun=mean) )

If you simply want to subset the pixels, as a point feature class, you can coerce your raster to points and use over to subset them. You can use writeOGR in rgdal, to write a shapefile.

pts <- as(r, "SpatialPointsDataFrame")
  proj4string(polys) <- proj4string(r)
z <- pts[!is.na(sp::over(pts, sp::geometry(polys))), ]
z@data <- data.frame(z@data, stats::na.omit(sp::over(pts, polys)))
z@data <- data.frame(z@data, coordinates(z)) 
head(z@data)

plot(z,pch=20,cex=0.5)