如何将UTM坐标转换为R中的纬度和经度

Question

我正在尝试运行物种分布模型，需要创建背景点来运行我的逻辑回归模型。我刚刚创建了 500 个随机点，但它们位于 UTM 坐标中，我需要纬度和经度。有没有办法在 R 中将它们转换为纬度和经度？如果是这样，你能与我分享代码吗？我对 R 还很陌生。谢谢！

Answer 1

如果您需要长/纬度，您可能应该使用该坐标参考系统生成随机点。但除此之外，你也可以做如下的事情。

我首先生成一组示例坐标 (

points

)

 set.seed(1)
 x <- runif(10, -10000, 10000)   
 y <- runif(10, -10000, 10000)   
 points <- cbind(x, y)

现在，假设您知道点的坐标参考系 (CRS)，您可以创建一个空间对象并分配该已知的 CRS。在我的示例中，这些点位于

UTM, zone 10, datum=WGS84

投影中。

library(terra)
v <- vect(points, crs="+proj=utm +zone=10 +datum=WGS84  +units=m")
v
# class       : SpatVector 
# geometry    : points 
# dimensions  : 10, 0  (geometries, attributes)
# extent      : -8764.275, 8893.505, -6468.865, 9838.122  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 

现在我们可以将它们转换为另一个 CRS，例如

longitude/latitude

y <- project(v, "+proj=longlat +datum=WGS84")
y
# class       : SpatVector 
# geometry    : points 
# dimensions  : 10, 0  (geometries, attributes)
# extent      : -127.5673, -127.4091, -0.05834327, 0.08873723  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
 

你可以像这样提取坐标

lonlat <- geom(y)[, c("x", "y")]
head(lonlat, 3)
#             x             y
#[1,] -127.5308 -0.0530354276
#[2,] -127.5117 -0.0583432750
#[3,] -127.4757  0.0337371933

你当然也可以做相反的事情

back <- project(y, "+proj=utm +zone=10 +datum=WGS84 +units=m")

可以使用

sf

包或旧的

sp

包完成同样的操作。使用

sp

，创建一个

SpatialPoints

对象并使用

spTransform

。

library(sp)
sputm <- SpatialPoints(points, proj4string=CRS("+proj=utm +zone=10 +datum=WGS84")) 
spgeo <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
lnlt <- coordinates(spgeo)
 

我在示例中使用了 UTM 区域 10。但请注意，UTM 有 60 个区域，您必须选择一个。每个覆盖 (360/60=) 6 度的条带。如果您的数据跨越大量经度或跨越 UTM 区域，则不应使用 UTM。对于 [-180, 180) 之间的经度，您可以像这样计算您需要的区域

utm_zone <- function(longitude) { trunc((180 + longitude) / 6) + 1 }

longs <- c(-122,-119, -118)
utm_zone(min(longs))
# [1] 10
utm_zone(max(longs))
# [1] 11

utm_zone(max(longs))

或者看看像这个这样的地图

与 UTM 的一个混淆点是，在南半球的位置使用“假北距”以避免负坐标。这是通过将 10,000,000 添加到 y 坐标来完成的，如下图所示，使用

+south

元素。

s <- vect(cbind(174, -44), crs="+proj=longlat +datum=WGS84")
geom(project(s, "+proj=utm +zone=59"))[, c("x", "y")]
#       x          y 
#740526.3 -4876249.1 

geom(project(s, "+proj=utm +south +zone=59"))[, c("x", "y")]
#       x         y 
#740526.3 5123750.9 

另请注意，我使用“PROJ4”符号来定义 CRS。如果使用的基准（基准是地球表面形状的模型）是 WGS84 或 NAD83，则效果很好。如果不是这种情况，您将需要使用 CRS 的“EPSG”代码或“WKT2”描述。