Chuyển đổi vĩ độ và kinh độ thành UTM

Question

Tôi đã tìm thấy một ví dụ khá đơn giản về cách thực hiện điều này nhưng tôi không thể làm cho nó hoạt động với tôi. Tôi khá mới đối với 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 

Nhưng đây là số ví dụ. Tôi có hàng ngàn tọa độ tôi phải chuyển đổi và tôi không thể tìm ra cách để đưa chúng từ bảng của tôi vào tập lệnh này

Tập dữ liệu của tôi có 3 cột, ID, X và Y. phương trình? Tôi đã bị kẹt trong tuần này

Answer 1

Trong câu hỏi của bạn, bạn không rõ liệu bạn đã đọc trong tập dữ liệu của mình thành một data.frame hay ma trận. Vì vậy, tôi cho rằng trong những điều sau đây bạn có dữ liệu bạn đặt trong một tập tin văn bản:

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

Cập nhật:

Các ý kiến ​​cho rằng vấn đề xuất phát từ áp dụng project() để data.frame. project() không hoạt động trên data.frames vì ​​nó kiểm tra is.numeric(). Do đó, bạn cần chuyển đổi dữ liệu thành ma trận như trong ví dụ trên. Nếu bạn muốn dính vào mã của bạn có sử dụng cbind() bạn phải làm như sau:

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

Sự khác biệt giữa dd["X"] và dd[,"X"] là sau này sẽ không quay trở lại một data.frame và kết quả là cbind() sẽ mang lại một ma trận và không phải là data.frame.

Answer 2

Để đảm bảo rằng siêu dữ liệu chiếu phù hợp ở mọi bước được liên kết với tọa độ, tôi khuyên bạn nên chuyển đổi điểm thành đối tượng SpatialPointsDataFrame càng sớm càng tốt.

Xem ?"SpatialPointsDataFrame-class" để biết thêm về cách chuyển đổi dữ liệu đơn giản.frames hoặc ma trận sang SpatialPointsDataFrame đối tượng.

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

Chuyển đổi vĩ độ và kinh độ điểm để 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

Dựa trên mã trên, tôi cũng đã thêm một phiên bản của khu vực và phát hiện bán cầu (mà giải quyết vấn đề chuyển đổi, như mô tả trong một số ý kiến) + viết tắt cho CRS chuỗi và chuyển đổi trở lại 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) 
} 

đâu CRS("+init=epsg:4326") là viết tắt cho CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0").

tìm tài liệu tham khảo khu UTM: http://www.igorexchange.com/node/927