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

r coordinates rgdal

42493 просмотра

5 ответа

71 Репутация автора

Я нашел довольно простой пример того, как это сделать, но не могу заставить его работать на меня. Я довольно новичок в 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. Как я могу преобразовать их, используя это уравнение? Я застрял на этом в течение нескольких недель

Автор: Colin Источник Размещён: 05.09.2013 03:09

Ответы (5)


4 плюса

25082 Репутация автора

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

# 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
# ...

Обновить:

Комментарии предполагают, что проблема возникает из-за применения 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.

Автор: f3lix Размещён: 05.09.2013 03:40

24 плюса

135212 Репутация автора

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

Подробнее ?"SpatialPointsDataFrame-class"о том, как преобразовать простые data.frames или матрицы в объекты 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 
Автор: Josh O'Brien Размещён: 05.09.2013 04:01

13 плюса

1823 Репутация автора

Преобразование точек широты и долготы в 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)
Автор: Stanislav Размещён: 13.05.2015 09:52

2 плюса

2334 Репутация автора

На основе приведенного выше кода я также добавил версию обнаружения зоны и полушария (которая решает проблемы преобразования, как описано в некоторых комментариях) + сокращение для 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

Автор: Teodor Ciuraru Размещён: 17.02.2018 05:16

0 плюса

1553 Репутация автора

Что касается примера, то зона UTM по умолчанию для данных координат составляет 50. Не рекомендуется проецировать координаты в дальние зоны. Вы можете проверить конвертацию с помощью инструмента NCAT от NOAA.

Код ниже использует sfпакет для преобразования.

library(sf)
library(tidyverse)

# Coordinate examples with expected UTM values
coord_sample <- data.frame(
  "Northing" = c(1105578.589, 5540547.370),
  "Easting" = c(609600.773, 643329.124),
  "Latitude" = c(10, 50),
  "Longitude" = c(118, 119))

#' Find UTM EPSG code from Latitude and Longitude coordinates (EPSG 4326 WGS84)
#' (vectorised)
#' Source: https://geocompr.robinlovelace.net/reproj-geo-data.html
#' Source: https://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point
LatLonToUTMEPSGCode <- function(lat, lon) {

  zone_number <- (floor((lon + 180) / 6) %% 60) + 1

  # Special zones for Norway
  cond_32 <- lat >= 56.0 & lat < 64.0 & lon >= 3.0 & lon < 12.0
  zone_number[cond_32] <- 32

  # Special zones for Svalbard
  cond_lat <- lat >= 72.0 & lat < 84.0

  cond_31 <- cond_lat & lon >= 0.0 & lon <  9.0
  zone_number[cond_31] <- 31

  cond_33 <- cond_lat & lon >= 9.0 & lon < 21.0
  zone_number[cond_33] <- 33

  cond_35 <- cond_lat & lon >= 21.0 & lon < 33.0
  zone_number[cond_35] <- 35

  cond_37 <- cond_lat & lon >= 33.0 & lon < 42.0
  zone_number[cond_37] <- 37

  # EPSG code
  utm <- zone_number
  utm[lat > 0] <- utm[lat > 0] + 32600
  utm[lat <= 0] <- utm[lat <= 0] + 32700

  return(utm)
}

sf_sample <- sf::st_as_sf(coord_sample, coords = c("Longitude", "Latitude"),
                          crs = 4326)

sf_sample %>%
  do(cbind(., st_coordinates(.))) %>%
  mutate(EPSG = LatLonToUTMEPSGCode(lat = Y, lon = X)) %>%
  group_by(EPSG) %>%
  do(cbind(as.data.frame(.) %>% select(Northing, Easting),
           st_coordinates(st_transform(., crs = head(.$EPSG, 1))))) %>%
  ungroup()


Вы можете проверить конверсию, сравнив с ожидаемыми значениями:

# A tibble: 2 x 5
   EPSG Northing Easting      X       Y
  <dbl>    <dbl>   <dbl>  <dbl>   <dbl>
1 32650  1105579  609601 609601 1105579
2 32650  5540547  643329 643329 5540547
Автор: Emer Размещён: 21.04.2019 10:19
Вопросы из категории :
32x32