将 WSG84 转换为 UTM，rgdal 替代方案

Question

我有许多现有脚本，其中我使用了

project()

包中的

rgdal

将 WSG84 纬度/经度坐标转换为 UTM。

rgdal

在更新后随机停止工作，间歇性地让我头疼，我总是能够通过将旧版本复制到我的包文件夹而不是从 CRAN 安装新版本来解决这个问题。现在，这也已停止工作，并且由于

rgdal

将在 2023 年退休，我正在尝试使用不同的投影包转换我的代码。这是我卡住的地方，因为我不太明白如何将它调整为其他包，如

terra

或

sf

.

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"))

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 

Answer 1

{rgdal} 包现在已经被取代，包维护者（已经退休）正在积极推动包用户考虑其他选择。因此，大胖警告和指向 {sf} 和 {terra} 的指针。

为了将 Long Lat 数据转换为 UTM zone 4N 坐标，我建议使用以下代码。

关键部分是：

• 从常规数据框转换为 {sf} 数据框；这将为您提供投影的基线（您从哪里投影？您的原始坐标是十进制还是米？甚至测量英尺？）
• 将

  sf::st_transform()

  应用到您的 sf 对象，从已知的开始转换为所需的格式

我发现使用 EPSG 代码（例如 4326 和 32604）而不是代码中使用的 proj4strings 来指定坐标参考系是最简单的。

library(sf)
library(dplyr)

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

# this is the baseline - locations in WGS84
locations_sf <- somelocations %>% 
  st_as_sf(coords = c("Long", "Lat"), crs = 4326)

# a visual check - do the locations look correct?
mapview::mapview(locations_sf)

# this is the action!
locations_projected <- locations_sf %>% 
  st_transform(crs = "EPSG:32604") # see https://epsg.io/32604 for more info...

# do the coordinates look as expected?
locations_projected
# Simple feature collection with 3 features and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 742693.4 ymin: 2302727 xmax: 746948.4 ymax: 2307439
# Projected CRS: WGS 84 / UTM zone 4N
#   id                 geometry
# 1  X POINT (746948.4 2302727)
# 2  Y POINT (745008.7 2305036)
# 3  Z POINT (742693.4 2307439)

Answer 2

这里是你如何用“terra”做到这一点。

您的示例数据

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

“原始”方法，与 rgdal 最相似：

library(terra)
project(as.matrix(somelocations[,c('Long', 'Lat')]), 
        "+proj=longlat", "+proj=utm +zone=4 +units=m")
#         [,1]    [,2]
#[1,] 746948.4 2302727
#[2,] 745008.7 2305036
#[3,] 742693.4 2307439

或者通过 SpatVector（这可能对其他原因有用，例如数据分析和映射）。

v <- vect(somelocations, geom=c("Long", "Lat"), crs="+proj=longlat")
p <- project(v, "+proj=utm +zone=4 +units=m")
crds(p)
#         [,1]    [,2]
#[1,] 746948.4 2302727
#[2,] 745008.7 2305036
#[3,] 742693.4 2307439