Extract cell values and coordinates from raster using a polygon in R

Question

I'm trying to extract cell values and coordinates from raster using a polygon in R. I have multiple raster layers with different types of values, all within the same extent. I have multiple polygons within these rasters data set areas.

I want to use the polygons to extract the raster cell values and coordinates within each cell.

I'm using this code for example:

library("raster")
library("rgdal")
library("shapefiles")
library("rgeos")
library("sp")

values_cps <- extract(CPs,cp_layer,cellnumbers=T)

cp_values <- as.data.frame(values_cps)

write.csv(cp_values,file="RCP95_cps.csv")

where: CPs is a raster layer, cp_layer is a polygon.

When I write the csv it gives me a column of cell ID numbers and a column of the values of each cell. I want to get the coordinates to come with it.

I've tried using coordinates() which I'm successfully using for the raster file, but not the extracted file.

Any ideas?

Answer 1

you can use the raster::xyFromCell function to get a matrix of coordinates from a vector of cell numbers, which you can cbind onto the extract (as the xyFromCell result is in exactly the same order as the vector provided)

I've used sf but of course this works with readOGR etc.;

require(raster)
require(sf)

# make 2 dummy rasters for a dummy stack
r <- raster(xmn=0,xmx=10,ymn=0,ymx=10,res=1)
r[] <- sample(1:4,100,replace=T)

r2 <- raster(xmn=0,xmx=10,ymn=0,ymx=10,res=1)
r2[] <- sample(8:10,100,replace=T)

s <- stack(r,r2)

# make a triangle simple feature
s2 <- rbind(c(2.2,2.2), c(3.3,8.5), c(5.3,0.8), c(2.2,2.2))
s2.sf <-st_sfc(st_polygon(list(s2)))
s2.pol = st_sf(ID = "tr", s2.sf)

# extract values into data frame
ex.df <- as.data.frame(extract(s,s2.pol,cellnumbers=T))

# create coordinate columns using xyFromCell
ex.df.coords <- cbind(ex.df, xyFromCell(s,ex.df[,1]))

Answer 2

You can get a table of cell values, coverage fractions, and center coordinates with the include_xy argument to exactextractr::exact_extract. Here's an example:

library(raster)
library(sf)
library(exactextractr)

# Pull municipal boundaries for Brazil
brazil <- st_as_sf(getData('GADM', country='BRA', level=2))[1:10, ]
    
# Pull gridded precipitation data
prec <- getData('worldclim', var='prec', res=10)[[12]]
    
prec_dfs <- exact_extract(prec, brazil, include_xy=TRUE, include_cols='GID_2')

This produces a list (prec_dfs) with one data frame for each feature in brazil. Each data frame looks like this:

       GID_2 value         x          y coverage_fraction
1  BRA.1.1_1   227 -67.08333  -9.750000        0.42837274
2  BRA.1.1_1   223 -66.91667  -9.750000        0.39436889
3  BRA.1.1_1   219 -66.75000  -9.750000        0.03865873
4  BRA.1.1_1   230 -67.08333  -9.916667        0.60549742
5  BRA.1.1_1   225 -66.91667  -9.916667        1.00000000

where value is the value of prec at the cell whose center coordinates are x and y, coverage_fraction of which is covered by the polygon.

These data frames can then be combined into a single table with rbind:

tbl <- do.call(rbind, prec_dfs)