将纬度和经度点转换为UTM

问题描述：

我发现了一个相当简单的例子，说明如何做到这一点，但我不能让它为我工作。我很新的R将纬度和经度点转换为UTM

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.我怎样才能使用它方程？我一直坚持这个星期

这将是很难帮助你，如果你不给我们一些数字是**不**工作（加上可能的格式的一些描述它们被存储）。另外，你是否知道所有经纬度坐标都落在单个UTM区域？ –

那么它没有那么多，数字不起作用，因为我没有把数字。我需要脚本从数千个坐标表中读取它。我只是不知道该怎么做。 dd Colin

什么样的表？（它存储在一个文本文件？csv？关系数据库？一个Excel文件？等等等）等等）听起来你真的在这一点上问如何读取数据到R.在这种情况下，尝试使用Google或在SO搜索栏中搜索以上权限，然后再次询问您是否无法弄清楚。 –

答

在你的问题，你不清楚你是否已经读入你的数据集到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。由于它检查is.numeric()，因此project()不适用于data.frames。因此，您需要将数据转换为矩阵，如上例所示。如果你要坚持你的代码，使用cbind()你必须做到以下几点：

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

dd["X"]和dd[,"X"]之间的区别是，后者将不会返回一个data.frame并因此cbind()将产生一个矩阵而不是一个data.frame。

谢谢。这看起来似乎产生了一种解脱！输出是以米为单位的？我现在的大多数坐标看起来像-9596387 -5670257，-9597204 -5666367等。这看起来是否正确？有些人提到了UTM区域。从谷歌我建议似乎分布在39升，39 K和38 K如果这是有道理的。脚本中是否必须考虑到这一点。正如我所说，我刚刚复制了一个我在网上找到的例子。对不起，如果我听起来像一个白痴，我只是想抓住这一点。再次感谢 – Colin

答

为了确保适当的投影元数据在与坐标相关的每个步骤，我建议尽快将点转换为SpatialPointsDataFrame对象。

有关如何将简单数据帧或矩阵转换为SpatialPointsDataFrame对象的更多信息，请参阅?"SpatialPointsDataFrame-class"。

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

spTransform（）调用中的变量“SP”是什么？如果我将其更改为“xy”，则不会得到如下所示的相同输出。 – stackoverflowuser2010

@ stackoverflowuser2010 - 糟糕，我的不好。那个'SP'从我编辑后的初始答案中遗留下来，它涉及一个'SpatialPoints'对象而不是'SpatialPointsDataFrame'。为了您将来的参考，您可以通过点击已编辑答案正下方的“编辑...之前”链接查看早期版本的答案。感谢您的支持。 –

答

转换纬度和经度点到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)

答

基于上述代码中，我还添加一个版本区和半球检测（即解决了转换的问题，如在一些注释中描述）的+速记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")是速记10。

寻找UTM区参考：http://www.igorexchange.com/node/927

将纬度和经度点转换为UTM

相关推荐