Extract Raster Pixels Values Using Vector Polygons in R

user_13013341 0
Viewed 1021

I have been struggling with this for hours. I have a shapefile (called "shp") containing 177 polygons i.e. 177 counties. This shapefile is overlaid on a raster. My raster (called "ras") is made of pixels having different pollution values.

Now I would like to extract all pixel values and their number of occurrences for each polygon.

This is exactly what the QGIS function "zonal histogram" is doing. But I would like to do the exact same thing in R.

I tried the extract() function and I managed to get a mean value per county, which is already a first step, but I would like to make a pixels distribution (histogram).

Could someone give me a hand ?

Many thanks,

Marie-Laure

raster gis r
2 Answers

Thanks a lot for your help. Next time I promise I will be careful and explain my issue more in details.

With your help I managed to find a solution. I also used this website : http://zevross.com/blog/2015/03/30/map-and-analyze-raster-data-in-r/

For information, first I had to uninstall the "tidyr" package because there was a conflict with the extract function.

In case it can help someone, here is the final code :

# Libraries loading
library(raster) 
library(rgdal)
library(sp)

# raster layer import
ras=raster("C:/*.tif")

# shapefile layer import
shp<-shapefile("C:/*.shp")

# Extract the values of the pixels raster per county
ext <- extract(ras, shp, method='simple')

# Function to tabulate pixel values by region & return a data frame
tabFunc                            <- function(indx, extracted, region, regname) {
  dat                              <- as.data.frame(table(extracted[[indx]]))
  dat$name                         <- region[[regname]][[indx]]
  return(dat)
}

# run through each county & compute a table of the number
# of raster cells by pixel value. ("CODE" is the county code) 
tabs <- lapply(seq(ext), tabFunc, ext, shp, "CODE")

# assemble into one data frame
df <- do.call(rbind, tabs)  

# to see the data frame in R
print(df)

# table export 
write.csv(df,"C:/*.csv", row.names = FALSE)
user_130133410

Here is a minimal, self-contained, reproducible example (almost literally from ?raster::extract, so not difficult to make)

library(raster)
r <- raster(ncol=36, nrow=18, vals=rep(1:9, 72))
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- spPolygons(cds1, cds2)

Now you can do

v <- extract(r, polys)
par(mfrow=c(1,2))
z <- lapply(v, hist)

Or more fancy

mains <- c("first", "second")
par(mfrow=c(1,2))
z <- lapply(1:length(v), function(i) hist(v[[i]], main=mains[i]))

Or do you want a barplot

z <- lapply(1:length(v), function(i) barplot(table(v[[i]]), main=mains[i]))
user_6352450
Related