基于坐标提取生物地理区域时R代码报错求助

我希望根据经纬度坐标提取各采样点对应的欧洲生物地理区域，下载对应shapefile后运行代码，却遇到identicalCRS(x, y) is not TRUE错误。

原加载shapefile代码 #

pg <- rgdal::readOGR(dsn = "BiogeoRegions2016.shp",
                     layer = "BiogeoRegions2016", p4s = "+init=epsg:3035") #load data
proj4string(pg) #Check projection
#plot(pg) #Plot

原示例点处理代码 #

# example points
pt <- data.frame(x = c(14,28), y = c(41, 59), id = 1:2) #example points
coordinates(pt) <- ~x+y #Transform to spatial coordinates
proj4string(pt) <- CRS("+init=epsg:4326") #projection
proj4string(pt) #projection

#  project to CRS of pg
pt <- spTransform(pt, CRS("+init=epsg:3035")) #Transform the point to the same projection
over(pt, pg) #This function should report the biogeographical region but I got the next error:

错误信息 #

Error in .local(x, y, returnList, fn, ...) :
identicalCRS(x, y) is not TRUE

解决方案 #

问题根源 #

手动指定的投影参数（p4s = "+init=epsg:3035"）和shapefile自带的投影信息可能存在细微格式差异，导致identicalCRS()判断不通过。shapefile本身已包含完整投影信息，无需手动指定。

修正后的sp体系代码 #

# 读取shapefile，不手动指定投影
pg <- rgdal::readOGR(dsn = "BiogeoRegions2016.shp", layer = "BiogeoRegions2016")
proj4string(pg) # 确认投影

# 处理示例点
pt <- data.frame(x = c(14,28), y = c(41, 59), id = 1:2)
coordinates(pt) <- ~x+y
proj4string(pt) <- CRS("+init=epsg:4326")

# 直接使用pg的投影参数进行转换，避免手动输入误差
pt <- spTransform(pt, proj4string(pg))
# 执行空间关联，获取对应区域信息
over(pt, pg)

更简洁的sf包替代方案（推荐） #

现代空间数据处理推荐使用sf包，代码更直观且兼容性更好：

library(sf)

# 读取生物地理区域shapefile
pg_sf <- st_read("BiogeoRegions2016.shp")

# 创建示例点的sf对象，指定WGS84坐标系（EPSG:4326）
pt_sf <- st_as_sf(data.frame(x = c(14,28), y = c(41, 59), id = 1:2),
                  coords = c("x", "y"), crs = 4326)

# 转换投影并执行空间连接，直接得到每个点对应的区域信息
st_join(st_transform(pt_sf, st_crs(pg_sf)), pg_sf)

内容的提问来源于stack exchange，提问作者Isma Soto Almena