替代已退役rgdal包，为移动数据添加区域列的方法

问题背景 #

我有一份移动数据集（坐标采用EPSG:5321投影），示例如下：

longitude latitude DateTimeRounded
459663.3  7181890      2007-09-10
459734.2  7181938      2007-09-11
459680.5  7181933      2007-09-12
459640.1  7181893      2007-09-13
459605.2  7181897      2007-09-14
459928.7  7182175      2007-09-15
459855.1  7182104      2007-09-16

同时拥有一份加拿大省份边界的Shapefile（.shp格式），需要给数据集新增region列，标记每个点所属的省份，不在省份范围内的点标记为NA，最终效果示例：

longitude latitude DateTimeRounded  region
459663.3  7181890      2007-09-10   Manitoba
459734.2  7181938      2007-09-11   Manitoba
459680.5  7181933      2007-09-12   Manitoba
459640.1  7181893      2007-09-13   NA
459605.2  7181897      2007-09-14   NA
459928.7  7182175      2007-09-15   NA
459855.1  7182104      2007-09-16   Ontario

原使用已退役的rgdal包实现该功能，现需用terra或sp包替代。

Shapefile分享建议 #

Shapefile并非单一文件，而是由.shp（几何数据）、.shx（索引）、.dbf（属性表）等多个关联文件组成。若需分享或展示，需将所有相关文件打包成压缩包（如ZIP），并说明文件的投影信息（比如示例中的EPSG:5321）。

解决方案 #

方法1：使用terra包（推荐，当前主流空间处理包） #

terra是rgdal/raster的替代包，API简洁高效：

# 加载包
library(terra)

# 读取加拿大省份Shapefile
setwd("CanadaMap/")
can_map <- vect("canada 5321.shp")  # 直接读取shp文件，自动识别关联文件

# 将移动数据转换为SpatVector对象
df_sp <- vect(df, geom = c("longitude", "latitude"), crs = "EPSG:5321")

# 若Shapefile和数据投影不一致，执行转换；示例中投影一致，可跳过
# df_sp <- project(df_sp, crs(can_map))

# 执行空间关联，获取每个点所属省份
region_info <- extract(can_map, df_sp, "NAME")

# 将结果合并到原数据集
df$region <- region_info$NAME

# 查看结果
head(df)

方法2：使用sp包（兼容原有代码逻辑） #

sp包仍可使用，用sf替代rgdal完成Shapefile读取：

# 加载包
library(sp)
library(sf)

# 读取加拿大省份Shapefile
setwd("CanadaMap/")
can_map_sf <- st_read("canada 5321.shp")
can_map <- as(can_map_sf, "Spatial")  # 转换为sp对象

# 将移动数据转换为SpatialPointsDataFrame
coordinates(df) <- ~longitude + latitude
proj4string(df) <- CRS("EPSG:5321")

# 若Shapefile和数据投影不一致，执行转换；示例中投影一致，可跳过
# df <- spTransform(df, proj4string(can_map))

# 执行空间关联
df$region <- over(df, can_map)$NAME

# 转换回普通数据框
df <- as.data.frame(df)

# 查看结果
head(df)

内容的提问来源于stack exchange，提问作者Cam