WGS84转UTM：rgdal替代方案及代码适配难题

原问题与报错 #

原有大量脚本依赖rgdal::project()完成WGS84经纬度到UTM坐标的转换，但rgdal已在2023年退役，更新后频繁出现运行故障，报错信息如下：

Error in rgdal::project(as.matrix(DASARs[, c("Long", "Lat")]), " +proj=utm +zone=4N ellps=WGS84 +units=m") :
target crs creation failed: Invalid value for an argument
In addition: Warning message:
PROJ support is provided by the sf and terra packages among others

原测试代码：

somelocations<-data.frame(id=c('X', "Y", "Z"),
                   Long=c(-156.6274,-156.6457,-156.6676),
                   Lat=c(20.8081,20.8292,20.8512))

require(rgdal)

rgdal::project(as.matrix(somelocations[,c('Long', 'Lat')]), " +proj=utm +zone=4N ellps=WGS84 +units=m"))

terra包替代方案 #

terra是rgdal的官方替代工具之一，支持无空间对象的直接坐标转换，兼容性更强：

# 安装并加载terra
install.packages("terra")
library(terra)

# 直接转换坐标矩阵
coords <- as.matrix(somelocations[, c("Long", "Lat")])
utm_coords <- project(coords, from = "EPSG:4326", to = "EPSG:32604")

# 合并转换结果到原数据框
somelocations <- cbind(somelocations, UTM_E = utm_coords[,1], UTM_N = utm_coords[,2])

若需要处理空间对象，可使用更完整的流程：

# 转为SpatVector空间对象
vec <- vect(somelocations, geom = c("Long", "Lat"), crs = "EPSG:4326")
# 投影转换
vec_utm <- project(vec, "EPSG:32604")
# 提取坐标
somelocations$UTM_E <- x(vec_utm)
somelocations$UTM_N <- y(vec_utm)

sf包替代方案 #

sf包更适合面向空间数据框的操作，语法直观：

# 安装并加载sf
install.packages("sf")
library(sf)

# 转为sf空间数据框
sf_locations <- st_as_sf(somelocations, coords = c("Long", "Lat"), crs = 4326)
# 投影转换
sf_utm <- st_transform(sf_locations, crs = 32604)
# 提取坐标到原数据框
somelocations$UTM_E <- st_coordinates(sf_utm)[,1]
somelocations$UTM_N <- st_coordinates(sf_utm)[,2]

关键注意事项 #

• 优先使用EPSG代码（如EPSG:4326代表WGS84，EPSG:32604代表UTM Zone 4N）替代手动编写PROJ字符串，避免格式错误
• terra和sf均支持最新PROJ标准，不会出现rgdal的兼容性问题
• 批量修改脚本时，可全局替换rgdal::project()调用为对应工具的转换逻辑

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