Extracting variables from grid cells using R

Question

I am trying to extract bioclimatic variables from individual grid cells. So i started loading my shapefile:

# I used the Admin 0 - Countries (1:10) dataset from Natural Earth 
shp <- sf::st_read(""ne_10m_admin_0_countries.shp"")
shape <- as_Spatial(shp)
crs <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
proj4string(shape) = crs

I used spsample() and HexPoints2SpatialPolygons() to create a hexagonal grid based on my shape object, and then I intersected the grid and the polygon.

size <- 0.5 #0.5 degrees
hex_points <- spsample(shape, type = "hexagonal", cellsize = size) 
hex_grid <- HexPoints2SpatialPolygons(hex_points, dx = size)
shape.grid <- gIntersection(shape, hex_grid, byid = T)

I load my points (646) and plotted then on my new shapefile: *Note that some of my points have the same coordinates.

coords <- read.table("clipboard", header=T)
coords$lat <- as.numeric(coords$lat)
coords$long <- as.numeric(coords$long)
coordinates(coords) <- ~long + lat
prj<-'+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'
coords <- SpatialPoints(coords, proj4string = CRS(prj))
proj4string(coords)<-proj4string(shape.grid)

Then I selected the grid cells which have points on (261): And here I am getting some problems.

shape.grid.containing.points <- shape.grid[coords, ]
bio10 <- raster(choose.files()) #[Bioclimatic variable][3] from WorldClim
bio10 <- extract(x=bio10, y=shape.grid.containing.points, fun=mean, cells = T)

Indeed I was able to extract variables from the grid cells i was aiming. However, the number of grid cells is smaller than the number of points, and I dont know which points were ommited because they where in the same grid or have the same coordinate.

How do I have a data frame with all of my points (646 lines) and the value of the variable I extract from the grid cell, instead of a data frame with only 261?

Answer 1

This might be a little bit confusing now, but let me provide an answer connected to your last question - How can i extract mean bioclimatic variables from unprojected grid cells? - using sf.

You're not getting the result desired, since data_sf contains 60 elements but grid_subset only has 11 after st_filter() because multiple hits are not considered:

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
library(rgbif)

shp <- st_read("ne_10m_admin_0_countries.shp")
#> Reading layer `ne_10m_admin_0_countries' from data source 
#>   `ne_10m_admin_0_countries.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 258 features and 168 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
#> Geodetic CRS:  WGS 84

grid <- st_make_grid(shp,
                     cellsize = 2.5,
                     crs = 4326,
                     square = TRUE) |> st_as_sf()

gbif_data <- occ_data(scientificName = 'Lestes sponsa',
                      hasCoordinate = TRUE, 
                      limit = 60)
gbif_data <- gbif_data$data

data_sf <- st_as_sf(gbif_data, 
                    coords = c("decimalLongitude", "decimalLatitude"), 
                    crs = st_crs(4326))

dim(data_sf)
#> [1]  60 102

grid_subset <- st_filter(grid, data_sf)

dim(grid_subset)
#> [1] 11  1

There are probably many approaches to solve this; here is one I tried instead of using st_filter() and which gave me plausible results from my point of view. st_intersects() checks for the topological relation (resp. geometric predicate) "intersects" between your inputs. The returned indices are used to query your grid cells. Your result consists of 60 grid cells using this, thereof 11 unique ones.

grid_subset <- grid[st_intersects(data_sf, grid) |> unlist(), ]

grid_subset 
#> Simple feature collection with 60 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -2.5 ymin: 45 xmax: 20 ymax: 62.5
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                                 x
#> 1  POLYGON ((10 57.5, 12.5 57....
#> 2  POLYGON ((10 57.5, 12.5 57....
#> 3  POLYGON ((-2.5 47.5, 8.5265...
#> 4  POLYGON ((10 45, 12.5 45, 1...
#> 5  POLYGON ((5 50, 7.5 50, 7.5...
#> 6  POLYGON ((5 50, 7.5 50, 7.5...
#> 7  POLYGON ((12.5 47.5, 15 47....
#> 8  POLYGON ((5 50, 7.5 50, 7.5...
#> 9  POLYGON ((5 50, 7.5 50, 7.5...
#> 10 POLYGON ((5 50, 7.5 50, 7.5...

You should still be able to use this with terra::extract() in your next step.