使用R中多个UTM区域的点将纬度和经度数据转换为UTM

Question

此问题是此线程中所问问题的扩展：Converting latitude and longitude points to UTM。我不确定是否应该在那里回复，但我现在在这里问。

我有一个包含生物记录的数据集，其中 lat/long 坐标来自多个 UTM 区域和各种属性（日期、记录器、数量等）。我想将这些转换为 UTM 坐标。当点都在同一个 UTM 区域时，我找到了将 lat/long 转换为 UTM 的方法，但是当它们都在不同区域时，我正在努力寻找有效的解决方案。我能找到的唯一问类似问题的问题是这个问题：Projecting long/lat with multiple UTM zones，建议 OP 做一些完全不同的事情。

我已经尝试了（第一个）链接线程中的答案，其中有一个用于查找 long/lat 坐标的 UTM 带的函数，以及一个用于将 long/lat 坐标转换为基于指定区域的UTM。

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
}

LongLatToUTM <-
  function(x, y, zone) {    
    xy <- data.frame(X = x, Y = y, zone = zone) #Makes a dataframe with a column x coordinates, y coordinates and UTM zone.
    coordinates(xy) <- c("X", "Y") #generates coordinates in that dataframe
    proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  #sets projection (it's currently standard lat long)
    res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))  #changes the projection to UTM, based on the zone in the column.
    return(res)
  }

在我看来，它们应该在 lat/long 坐标列上工作正常，但第一个只是 returns 每行第一个坐标的 UTM 区域，第二个转换为第一个坐标的 UTM 区域（因此只有这个坐标和恰好位于同一区域的任何后续坐标是正确的）。显然我对这些工具还不够了解...

有人对如何阻止这种情况有任何建议吗？我觉得这可能更多的是关于语法和数据帧的问题，而不是处理空间数据，对这些函数的编写方式进行简单调整就可以解决问题。

此处的一些示例数据包含纬度 (x) 和经度 (y) 坐标、'correct' UTM 等价物（x_correct 和 y_correct）、它们正确的 UTM 区域，以及他们所在的国家。


test_coordinates <- data.frame(x = c(13.637398, -3.58627, -5.178889), y = c(41.30736, 40.72913, 40.17528), x_correct = c(385936, 450492, 314480), y_correct = c(4573773, 4508854, 4449488 ),  zone = c(33, 30, 30), key = c(1, 2, 3), country = c("italy", "spain", "spain"))

我的解决方案是使用 mapply 运行数据框上的 UTM 区域查找器，并编写一个新版本的转换器，依次获取每个点，运行也使用 mapply .

LongLatToUTM2 <-
  function(x, y, zone) {
    pt <- st_sfc(st_point(c(x, y), dim = XY)) #Creates a point out of x and y coordinates
    st_crs(pt) <- "+proj=longlat +datum=WGS84" #Sets the CRS as standard lat long for this point
    pt2<- st_transform(pt, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep=''))) #Converts this to UTM for the zone specified in the function
    return(pt2) #returns the new UTM point
  }

mapply(find_UTM_zone, test_coordinates$x, test_coordinates$y) -> zones2 #Creates a vector of the correct UTM zones

mapply(LongLatToUTM2, test_coordinates$x, test_coordinates$y, test_coordinates$zone) -> converted_points #Makes a list of the correct UTM points

data.frame(t(sapply(converted_points,c)))  #Makes this list into a dataframe

问题是速度很慢（我有 60,000 条记录要转换......）而且你最终得到一个 UTM 坐标数据框，你无法真正匹配原始坐标 -纵坐标（我不确定您是否真的可以将属性添加到 R 中的单个点？）。理想情况下，我希望能够保留具有转换后坐标的属性。

编辑：只是为了澄清我所追求的——我想改进这些功能，以便它们在数据帧上工作，而不是单一坐标，即输入经度和纬度数据帧，我可以添加一个包含其 UTM 区域的列，然后根据这些区域添加具有正确 UTM 坐标的另一列。

我意识到使用 dplyr 方法获取 UTM 区域相当容易（如果我现在忽略 Svalbard/Norway 坐标）：

test_coordinates %>% 
  mutate(zone2 = (floor((x + 180)/6) %% 60) + 1)

但是第二部分，如何轻松编写一个使用这些区域获取 UTM 坐标的函数对我来说仍然是个谜。

Answer 1

更新

这是一个使用 dplyr 和 spTransfrom

的更快更优雅的解决方法

增强数据（60k+ 行）：

test_coordinates <- data.frame(x = c(13.637398, -3.58627, -5.178889), y = c(41.30736, 40.72913, 40.17528), x_correct = c(385936, 450492, 314480), y_correct = c(4573773, 4508854, 4449488 ),  zone = c(33, 30, 30), key = c(1, 2, 3), country = c("italy", "spain", "spain"))
test_coordinates = rbind(test_coordinates, test_coordinates[rep(1,60*1000),]) # simulate big data

library(dplyr)
library(sp)

get_utm <- function(x, y, zone, loc){
  points = SpatialPoints(cbind(x, y), proj4string = CRS("+proj=longlat +datum=WGS84"))
  points_utm = spTransform(points, CRS(paste0("+proj=utm +zone=",zone[1]," +ellps=WGS84")))
  if (loc == "x") {
    return(coordinates(points_utm)[,1])
  } else if (loc == "y") {
    return(coordinates(points_utm)[,2])
  }
}

test_coordinates %<>% 
  mutate(zone2 = (floor((x + 180)/6) %% 60) + 1, keep = "all"
  ) %>% 
  group_by(zone2) %>% 
  mutate(utm_x = get_utm(x, y, zone2, loc = "x"),
         utm_y = get_utm(x, y, zone2, loc = "y"))

输出（仅 5 行）

test_coordinates


# A tibble: 603 × 10
# Groups:   zone2 [2]
       x     y x_correct y_correct  zone   key country zone2   utm_x    utm_y
   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl> <chr>   <dbl>   <dbl>    <dbl>
 1 13.6   41.3    385936   4573773    33     1 italy      33 385936. 4573773.
 2 -3.59  40.7    450492   4508854    30     2 spain      30 450492. 4508854.
 3 -5.18  40.2    314480   4449488    30     3 spain      30 314480. 4449488.
 4 13.6   41.3    385936   4573773    33     1 italy      33 385936. 4573773.
 5 13.6   41.3    385936   4573773    33     1 italy      33 385936. 4573773.

原始回复

替换：

data.frame(t(sapply(converted_points,c)))  #Makes this list into a dataframe

有：

test_coordinates$utm_x <- unlist(converted_points)[c(T,F)]
test_coordinates$utm_y <- unlist(converted_points)[c(F,T)]

          x        y x_correct y_correct zone key country    utm_x   utm_y
1 13.637398 41.30736    385936   4573773   33   1   italy 385935.7 4573773
2 -3.586270 40.72913    450492   4508854   30   2   spain 450492.4 4508854
3 -5.178889 40.17528    314480   4449488   30   3   spain 314479.5 4449488

使用R中多个UTM区域的点将纬度和经度数据转换为UTM

Converting latitude and longitude data to UTM with points from multiple UTM zones in R

r

spatial

coordinates

map-projections

sp