Как преобразовать координаты UTM в широту и долготу в R

Если вам нужна долгота/широта, вам, вероятно, следует генерировать случайные точки, используя эту систему отсчета координат. Но в противном случае вы можете сделать что-то вроде ниже.

Сначала я генерирую примерный набор координат (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. Если это не так, вам нужно будет использовать код EPSG или описание WKT2 вашей CRS.