Преобразование точек широты и долготы в UTM

Question

Я нашел довольно простой пример того, как это сделать, но я не могу заставить его работать для меня. Я довольно новичок к R

library(rgdal) 
xy <- cbind(c(118, 119), c(10, 50)) 
project(xy, "+proj=utm +zone=51 ellps=WGS84") 
      [,1] [,2] 
[1,] -48636.65 1109577 
[2,] 213372.05 5546301 

Но это номера примеров. У меня есть тысячи координат, которые мне нужно преобразовать, и я не могу понять, как их получить из таблицы в этот скрипт.

В моем наборе данных есть 3 столбца, ID, X и Y. Как я могу их преобразовать с помощью этого уравнение? Я застрял на этом в течение недель

Answer 1

В вашем вопросе вам не ясно, читаете ли вы данные в своем наборе данных в data.frame или matrix. Таким образом, я предполагаю, что в следующем у вас есть ваш набор данных в текстовом файле:

# read in data 
dataset = read.table("dataset.txt", header=T) 

# ... or use example data 
dataset = read.table(text="ID X Y 
1 118 10 
2 119 50 
3 100 12 
4 101 12", header=T, sep=" ") 

# create a matrix from columns X & Y and use project as in the question 
project(as.matrix(dataset[,c("X","Y")]), "+proj=utm +zone=51 ellps=WGS84") 
#    [,1] [,2] 
# [1,] -48636.65 1109577 
# [2,] 213372.05 5546301 
# ... 

Update:

Комментарии свидетельствуют о том, что проблема исходит от применения project() к data.frame. project() не работает на data.frames, так как он проверяет на is.numeric(). Поэтому вам нужно преобразовать данные в матрицу, как в моем примере выше. Если вы хотите придерживаться своего кода, который использует cbind() вы должны сделать следующее:

X <- dd[,"X"] 
Y <- dd[,"Y"] 
xy <- cbind(X,Y) 

Разница между dd["X"] и dd[,"X"] в том, что последний не будет возвращать data.frame и как следствие cbind() даст матрицу а не data.frame.

Answer 2

Для обеспечения того, чтобы соответствующие метаданные проекции находились на каждом шаге, связанном с координатами, я бы предложил как можно скорее преобразовать точки в объект SpatialPointsDataFrame.

См. ?"SpatialPointsDataFrame-class" подробнее о том, как конвертировать простые данные. Кадры или матрицы в SpatialPointsDataFrame объектов.

library(sp) 
library(rgdal) 

xy <- data.frame(ID = 1:2, X = c(118, 119), Y = c(10, 50)) 
coordinates(xy) <- c("X", "Y") 
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example 

res <- spTransform(xy, CRS("+proj=utm +zone=51 ellps=WGS84")) 
res 
#   coordinates ID 
# 1 (-48636.65, 1109577) 1 
# 2 (213372, 5546301) 2 

## For a SpatialPoints object rather than a SpatialPointsDataFrame, just do: 
as(res, "SpatialPoints") 
# SpatialPoints: 
#    x  y 
# [1,] -48636.65 1109577 
# [2,] 213372.05 5546301 
# Coordinate Reference System (CRS) arguments: +proj=utm +zone=51 
# +ellps=WGS84 

Answer 3

Преобразование широты и долготы точки для UTM

library(sp) 
library(rgdal) 

#Function 
LongLatToUTM<-function(x,y,zone){ 
xy <- data.frame(ID = 1:length(x), X = x, Y = y) 
coordinates(xy) <- c("X", "Y") 
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example 
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep=''))) 
return(as.data.frame(res)) 
} 

# Example 
x<-c(-94.99729,-94.99726,-94.99457,-94.99458,-94.99729) 
y<-c(29.17112, 29.17107, 29.17273, 29.17278, 29.17112) 
LongLatToUTM(x,y,15) 

Answer 4

На основании приведенного выше кода, я также добавил версию зоны и обнаружения полушария (который решает проблемы преобразования, как описано в некоторых комментариев) + стенография для CRS строки и преобразования обратно в WSG86:

library(dplyr) 
library(sp) 
library(rgdal) 
library(tibble) 

find_UTM_zone <- function(longitude, latitude) { 

# Special zones for Svalbard and Norway 
    if (latitude >= 72.0 && latitude < 84.0) 
     if (longitude >= 0.0 && longitude < 9.0) 
       return(31); 
    if (longitude >= 9.0 && longitude < 21.0) 
      return(33) 
    if (longitude >= 21.0 && longitude < 33.0) 
      return(35) 
    if (longitude >= 33.0 && longitude < 42.0) 
      return(37) 

    (floor((longitude + 180)/6) %% 60) + 1 
} 


find_UTM_hemisphere <- function(latitude) { 

    ifelse(latitude > 0, "north", "south") 
} 

# returns a DF containing the UTM values, the zone and the hemisphere 
longlat_to_UTM <- function(long, lat, units = 'm') { 

    df <- data.frame(
     id = seq_along(long), 
     x = long, 
     y = lat 
    ) 
    sp::coordinates(df) <- c("x", "y") 

    hemisphere <- find_UTM_hemisphere(lat) 
    zone <- find_UTM_zone(long, lat) 

    sp::proj4string(df) <- sp::CRS("+init=epsg:4326") 
    CRSstring <- paste0(
     "+proj=utm +zone=", zone, 
     " +ellps=WGS84", 
     " +", hemisphere, 
     " +units=", units) 
    if (dplyr::n_distinct(CRSstring) > 1L) 
     stop("multiple zone/hemisphere detected") 

    res <- sp::spTransform(df, sp::CRS(CRSstring[1L])) %>% 
     tibble::as_data_frame() %>% 
     dplyr::mutate(
      zone = zone, 
      hemisphere = hemisphere 
     ) 

    res 
} 

UTM_to_longlat <- function(utm_df, zone, hemisphere) { 

    CRSstring <- paste0("+proj=utm +zone=", zone, " +", hemisphere) 
    utmcoor <- sp::SpatialPoints(utm_df, proj4string = sp::CRS(CRSstring)) 
    longlatcoor <- sp::spTransform(utmcoor, sp::CRS("+init=epsg:4326")) 
    tibble::as_data_frame(longlatcoor) 
} 

Где CRS("+init=epsg:4326") представляет собой сокращенную CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0").

Поиск ссылки UTM зона: http://www.igorexchange.com/node/927