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:

DEM layer and plot shape overlay

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?


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.


    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(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))
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)) 

  • Jeffrey thanks for your comment. The v <- extract(r, polys) does not return any XY information of the pixels, which is what I want.
    – SinghD
    May 18 '17 at 16:04
  • You are not at all clear on this! It sounds like you just want the pixels as points intersecting your polygons, correct? I added to my answer to include this. May 18 '17 at 17:43
  • Let's say you had your polys IDs were 1001:1002. Now your 'v' list should include the 1001 index on pixels under that polygon 1001. With your original answer there is still no way to associate the polys real plot IDs to the points in raster.
    – SinghD
    May 18 '17 at 18:00
  • This is not the case, look at the list object. Each element in the list represents a polygon, with associated row.name, and the values are a vector of the raster values intersecting the polygon. If you use the "cellnumbers=TRUE" argument in raster::extract, then the explicit cell indices are returned. The list is also ordered the same as the polygons so, one can perform a one-to-one join of lapply results. May 18 '17 at 18:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.