Is there a way to get coordinates for raster values extracted from a polygon in R?

Question

I have a shapefile of polygons that I want to use to extract raster values into a data frame. So I do that in the following code.

shp <- sf:st_read('example.shp')
r  <- raster::raster('example.tif')

extract <- raster::extract(r, shp, df=TRUE)

This gives me a data frame of two columns: the numeric ID for each polygon and the associated extracted raster value. Now I would like to add the x, y coordinates for each extracted raster value. I have seen this done for point shapefiles but I am not sure how to apply it for polygon shapefile geometry.

I have just seen a very similar answer on stackoverflow.com/questions/48021657/… — dieghernan

dieghernan dieghernan · Accepted Answer · 2021-02-19T15:10:26

Try this, I combined raster::extract(..., cellnumbers = TRUE) and raster::xyFromCell:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
#> Loading required package: sp
library(giscoR)

# Dummy data
shp <-
  gisco_get_nuts(country = "Netherlands",
                 nuts_level = 3,
                 epsg = 4326)[, 1]
r <- raster::getData("alt", country = "Netherlands")
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
#> definition
extract <- raster::extract(r, shp, df = TRUE, cellnumbers = TRUE)
# Order (for checking purposes)
extract <- extract[order(extract$cell),]


# Extract coordinates
xy <- xyFromCell(r, cell = extract$cell, spatial = FALSE)

# Convert to df and add cellnumber
xy <- as.data.frame(xy)
xy$cell <- extract$cell

# Merge two data frames
extract_end <- merge(extract, xy)
extract_end <- extract_end[order(extract_end$cell),]

# This is what you are looking for: extract_end is a data frame
# has values and x and y are coordinates
head(extract_end)
#>   cell ID NLD_msk_alt        x        y
#> 1 2319  3          NA 6.620833 53.56250
#> 2 2796  3          NA 6.595833 53.55417
#> 3 2797  3          NA 6.604167 53.55417
#> 4 2798  3          NA 6.612500 53.55417
#> 5 2799  3          NA 6.620833 53.55417
#> 6 2800  3          NA 6.629167 53.55417

# Checks - can be ommited

nrow(extract) == nrow(extract_end)
#> [1] TRUE

# Check NAs in coordinates
nrow(extract_end[is.na(extract_end$x), ])
#> [1] 0
nrow(extract_end[is.na(extract_end$y), ])
#> [1] 0

# Convert to sf for checks

sfobj <- st_as_sf(extract_end, coords = c("x", "y"))
sfobj <- st_set_crs(sfobj, st_crs(4326))

# Plot as sf points
par(mar = c(3, 3, 3, 3))
plot(
  sfobj[, 3],
  axes = TRUE,
  main = "sf points",
  key.pos = 4,
  breaks = "equal",
  nbreaks = 100,
  pal = rev(terrain.colors(100))
)

# Compare with raster plot
par(mar = c(3, 3, 3, 3))
plot(r, main = "Raster")

sessionInfo()
#> R version 4.0.3 (2020-10-10)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19041)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
#> [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] giscoR_0.2.4-9000 raster_3.4-5      sp_1.4-5          sf_0.9-7         
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.6         pillar_1.4.7       compiler_4.0.3     highr_0.8         
#>  [5] class_7.3-18       tools_4.0.3        digest_0.6.27      evaluate_0.14     
#>  [9] lifecycle_1.0.0    tibble_3.0.6       lattice_0.20-41    pkgconfig_2.0.3   
#> [13] rlang_0.4.10       reprex_1.0.0       DBI_1.1.1          rgdal_1.5-23      
#> [17] yaml_2.2.1         xfun_0.21          e1071_1.7-4        styler_1.3.2      
#> [21] stringr_1.4.0      dplyr_1.0.4        knitr_1.31         generics_0.1.0    
#> [25] fs_1.5.0           vctrs_0.3.6        classInt_0.4-3     grid_4.0.3        
#> [29] tidyselect_1.1.0   glue_1.4.2         R6_2.5.0           rmarkdown_2.6     
#> [33] purrr_0.3.4        magrittr_2.0.1     codetools_0.2-18   backports_1.2.1   
#> [37] ellipsis_0.3.1     htmltools_0.5.1.1  units_0.6-7        assertthat_0.2.1  
#> [41] countrycode_1.2.0  KernSmooth_2.23-18 stringi_1.5.3      crayon_1.4.1

^{Created on 2021-02-19 by the reprex package (v1.0.0)}

Is there a way to get coordinates for raster values extracted from a polygon in R?

1 Answers