将纬度和经度点转换为UTM

I found a fairly simple example of how to do this but I cant get it to work for me. I'm pretty new to R

我找到了一个相当简单的例子来说明如何做到这一点，但是我不能让它为我工作。我对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

But this is with example numbers. I have thousands of coordinates I have to transform and I cant figure out how to get them from my table to into this script

但这是例子。我有成千上万的坐标，我必须变换，我不知道如何把它们从我的表中移到这个脚本中。

My data set has 3 columns, ID, X, and Y. How can I transform them using this equation? I've been stuck on this for weeks

我的数据集有3列，ID X y，我如何用这个方程来变换它们?我已经被这个问题困扰了好几个星期了

4 个解决方案

#1

To ensure that appropriate projection metadata are at every step associated with the coordinates, I'd suggest converting the points to a SpatialPointsDataFrame object as soon as possible.

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

See ?"SpatialPointsDataFrame-class" for more on how to convert simple data.frames or matrices to SpatialPointsDataFrame objects.

更多关于如何将简单数据、框架或矩阵转换为空间点数据aframe对象的信息，请参见“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 

#2

Converting latitude and longitude points to UTM

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

#3

In your question you are not clear whether you already read in your data set into a data.frame or matrix. So I assume in the following you have your data set in a text file:

在您的问题中，您不清楚您是否已经在数据集中读取到数据集或矩阵。因此，我假设您的数据集在以下文本文件中:

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

更新:

The comments suggest that the problem comes from applying project() to data.frame. project() does not work on data.frames since it checks for is.numeric(). Therefore, you need to convert data to a matrix as in my example above. If you want to stick to your code that uses cbind() you have to do the following:

注释表明问题来自于将project()应用到data.frame。project()对data.frame不起作用，因为它检查is.numeric()。因此，您需要将数据转换为上述示例中的一个矩阵。如果您想坚持使用cbind()的代码，您必须执行以下操作:

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

The difference between dd["X"] and dd[,"X"] is that latter will not return a data.frame and as a consequence cbind() will yield a matrix and not a data.frame.

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

#4

Based on the code above, I also added a version of zone and hemisphere detection (that solves conversion problems, as described in some comments) + shorthand for CRS string and conversion back to WSG86:

基于上面的代码，我还添加了一个版本的区域和半球检测(解决了转换问题，如某些注释所述)+ 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)
}

Where CRS("+init=epsg:4326") is shorthand for CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0").

其中CRS(“+init=epsg:4326”)是CRS的简写(“+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0”)。

Finding UTM zone reference: http://www.igorexchange.com/node/927

找到UTM区域参考:http://www.igorexchange.com/node/927。