我正在尝试运行物种分布模型,并且需要创建背景点来运行我的逻辑回归模型。我刚刚创建了 500 个随机点,但它们位于 UTM 坐标中,我需要经纬度。有没有办法将它们转换为 R 中的纬度和经度?如果是这样,你能把代码分享给我吗?我对R相当陌生。谢谢!
1 回答
如果您需要 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”描述。