Extract Raster Pixels Values Using Vector Polygons in R

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

Sounds challenging. I'm not sure that I can help, but there may be people in the community that can, even though they don't have a magic snippet in their utils folder. But that's unlikely without a reprex, since few will try to reverse-engineer a problem.

Please see the FAQ: What's a reproducible example (`reprex`) and how do I do one? Using a reprex, complete with representative data will attract quicker and more answers.

here's a brief post I wrote on the task, it uses sf and stars for the spatial polygon and raster data, respectively.

Hope if helps!

Glad you fixed it. Please mark the solution for the benefit of those to follow. No false modesty allowed.

For the data part of it: Only enough data to reproduce the issue is needed. It doesn't have to be your data, one of the built in datasets with comparable structure will do. It doesn't have to be all of the data, just enough to show the same behavior. And, as an added bonus.

dput(your_data

Will return a formal definition of the data object, which can be cut and pasted into the code reprex with

dput(your_data

followed by paste.

Here's a function with a slightly less verbose output:

require(clipr)
#> Loading required package: clipr
#> Welcome to clipr. See ?write_clip for advisories on writing to the clipboard in R.
require(magrittr)
#> Loading required package: magrittr
require(stringr)
#> Loading required package: stringr

specimen <- function(x)
  deparse(x) %>%
  str_c(collapse = '') %>%
  str_replace_all('\\s+', ' ') %>%
  str_replace_all('\\s*([^,\\()]+ =) (c\\()', '\n  \\1\n    \\2')  %>%
  str_replace_all('(,) (class =)', '\\1\n  \\2') %>%
  write_clip(allow_non_interactive = TRUE) 

(Credit to the anonymous author from whom I ripped off at least the major pieces. Mea culpa.)

Wah, I have plenty of work to learn how to use all those packages !

Thanks anyway

1 Like

Prepare yourself for a lifelong journey! Come on back with other questions as needed.

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

Hi technocrat,

Sorry about this. Next time I will use a reprex.

Is there a way we can attach a csv document here ?
Or is it better to use the tibble function ?

Many thanks for your reply

Meanwhile I found a solution to my issue.
First I had to uninstall the "tidyr" package because there was a conflict with the extract() function.
It took me ages to figure this out.
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)

Yes, thanks @Luis !

I forgot to mention that I used your website and also this one:

By the way, is there a way to find the data from Human Footprint index ? I did not manage to find it and stopped there !

The work you did is just amazing ! Many thanks

the human footprint dataset is on dryad :+1:
https://datadryad.org/stash/dataset/doi:10.5061/dryad.052q5

1 Like