如果你想从网格单元中提取平均值，你需要使用网格单元多边形而不是点坐标。你可以简单地选择覆盖点的多边形，然后提取这些多边形的平均栅格值，而不是使用“over”。

shape.grid.containing.points <- shape.grid[coords, ]

plot(shape.grid.containing.points)

bioclim_data <- extract(x=stackrasters, y=shape.grid.containing.points, fun=mean)

还需要注意的是GBIF数据通常需要清理，并且您不应该将CRS分配给已经有CRS的空间对象，如GADMMap。最终您需要迁移您的空间代码，例如迁移到'terra'，因为'sp'将被弃用。

请原谅我只使用类似的数据集来演示工作流（gadm_410约为1.4 GB，wc2.1_30s_bio约为9.7 GB）。我也会尽可能地坚持使用sf和terra：

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

# I used the Admin 0 - Countries (1:10) dataset from Natural Earth 
shp <- sf::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

# make hexagonal grid with res = 2.5° in WGS 84
grid <- sf::st_make_grid(shp,
                         cellsize = 2.5,
                         crs = 4326,
                         square = FALSE) |> sf::st_as_sf()

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

# create a simple features object from your data
data_sf <- sf::st_as_sf(gbif_data, 
                        coords = c("decimalLongitude", "decimalLatitude"), 
                        crs = sf::st_crs(4326))

# select objects from grid (= cells) containing points from data_sf (= locations)
grid_subset <- sf::st_filter(grid, data_sf)

几乎完成后，您只需导入栅格数据并使用terra::extract()来获得所需的值：

library(terra)

# I used wc2.1_30s_prec from WorldClim, read using `rast()`
files <- list.files(pattern = "*.tif")
prec <- terra::rast(files)

# note that the resulting SpatRast object has 12 layers
prec
#> class       : SpatRaster 
#> dimensions  : 21600, 43200, 12  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> sources     : wc2.1_30s_prec_01.tif  
#>               wc2.1_30s_prec_02.tif  
#>               wc2.1_30s_prec_03.tif  
#>               ... and 9 more source(s)
#> names       : wc2.1~ec_01, wc2.1~ec_02, wc2.1~ec_03, wc2.1~ec_04, wc2.1~ec_05, wc2.1~ec_06, ... 
#> min values  :           0,           0,           0,           0,           0,           0, ... 
#> max values  :         973,        1309,        1145,        1049,        2081,        2226, ... 

# extract values from prec by polygons from grid_subset using mean as aggregate
results <- terra::extract(prec, terra::vect(grid_subset), fun = mean)

# prec has 12 layers, grid_subset consists of 60 polygons, 
# c.f. dimensions below (+ID column containing an identifier of the related polygon)
str(results)
#> 'data.frame':    60 obs. of  13 variables:
#>  $ ID               : num  1 2 3 4 5 6 7 8 9 10 ...
#>  $ wc2.1_30s_prec_01: num  74 70 70 70 70 70 70 67 67 67 ...
#>  $ wc2.1_30s_prec_02: num  65 51 51 51 51 51 51 52 52 52 ...
#>  $ wc2.1_30s_prec_03: num  52 67 67 67 67 67 67 67 67 67 ...
#>  $ wc2.1_30s_prec_04: num  51 46 46 46 46 46 46 45 45 45 ...
#>  $ wc2.1_30s_prec_05: num  61 62 62 62 62 62 62 61 61 61 ...
#>  $ wc2.1_30s_prec_06: num  46 70 70 70 70 70 70 70 70 70 ...
#>  $ wc2.1_30s_prec_07: num  42 70 70 70 70 70 70 66 66 66 ...
#>  $ wc2.1_30s_prec_08: num  43 60 60 60 60 60 60 57 57 57 ...
#>  $ wc2.1_30s_prec_09: num  62 72 72 72 72 72 72 67 67 67 ...
#>  $ wc2.1_30s_prec_10: num  68 71 71 71 71 71 71 66 66 66 ...
#>  $ wc2.1_30s_prec_11: num  75 79 79 79 79 79 79 72 72 72 ...
#>  $ wc2.1_30s_prec_12: num  77 77 77 77 77 77 77 73 73 73 ...

对不起，我没有真正回答你的问题，而冒昧地投影了你的网格单元格。-）