Collectives™ on Stack Overflow

Find centralized, trusted content and collaborate around the technologies you use most.

Teams

Q&A for work

Connect and share knowledge within a single location that is structured and easy to search.

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.
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)
                Differences on plotting are because sf are points and raster are areas. You can check in the upper right corner, but values should be the same
– dieghernan
                Feb 19, 2021 at 15:14
        
Thanks for contributing an answer to Stack Overflow!
Please be sure to answer the question. Provide details and share your research!
But avoid …
Asking for help, clarification, or responding to other answers.
Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.