10

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

4

1 回答 1

24

如果您需要 long/lat，您可能应该使用该坐标参考系统生成随机点。但除此之外，您可以执行以下操作。

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

 library(terra)
 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(rgdal)
sputm <- SpatialPoints(points, proj4string=CRS("+proj=utm +zone=10 +datum=WGS84")) 
spgeo <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
lnlt <- coordinates(spgeo)

我在示例中使用了 UTM 区域 10。但请注意，有 60 个 UTM 区域，您必须选择一个。每个覆盖 (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”描述。

于 2015-05-03T19:47:50.047 回答