如何通过循环将分属不同UTM区的经纬度批量转换为UTM坐标

问题描述 #

我有一个包含纬度（lat）、经度（long）的数据框，需要将这些经纬度转换为UTM坐标，且这些经纬度分属不同的UTM投影带。数据框中已包含记录对应正确UTM带的utm_zone列，目前已有可转换单个UTM带经纬度的脚本，请问如何编写循环，依据各条数据对应的UTM带完成批量转换？

数据结构 #

data <- structure(list(site_code = c("pa1", "pa2", "pa5", "pa7", "pa4", "pa9", "pa11", "pa6", "pa12", "pa10", "pa3", "pa8", "tn1", "tn2", "tn5", "tn4", "tn3", "vt1", "vt2", "vt8", "vt7", "vt4", "vt3", "vt5", "vt6", "la2", "la1", "la4", "la3", "la5", "la6", "la7", "nm19", "nm10", "nm8", "nm4", "nm11", "nm1", "nm2", "nm15", "nm13", "nm14", "nm12", "nm6", "nm7", "nm3", "nm20", "nm5", "nm16", "nm18", "nm17", "nm36", "nm35", "nm9", "la8", "nm21", "nm22", "nm23", "nm24", "nm25", "nm26", "nm27", "nm28", "nm29", "nm30", "nm31", "nm32", "nm33", "nm34"), region = c("pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "pennsylvania", "tennessee", "tennessee", "tennessee", "tennessee", "tennessee", "vermont", "vermont", "vermont", "vermont", "vermont", "vermont", "vermont", "vermont", "louisiana", "louisiana", "louisiana", "louisiana", "louisiana", "louisiana", "louisiana", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "louisiana", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico", "new mexico"), site = c("wood lab pond", "rv pond", "tryon-weber", "beaver", "cow pit", "david's pond", "church", "tuttle pond", "black jack", "phelps pond", "hatchery", "vorisek pond", "newt", "pine", "sink", "deer", "west", "camp johnson pond 1", "camp johnson pond 2", "natural area close to vt01 and vt02", "adjacent to sand bar state park", "camp ethan allen pond", "north beach pond", "shelbourne bay", "shelbourne pond", "solomon lane", "kurthwood", "horse head permanent", "kisatchie bayou", "horse head ephemeral", "la6", "la7", "leasburg dam state park, bridge", "garcia well", "johnson tank", "davie's playa- doug burkett", "cuchillo", "red tank", "red tank 2", "bupu0039", "roadside ditch bude2026", "speamulti2237", "2330scacou", "rouse tank", "avilas", "rhodes spring", "davie's playa", "circle 7", "range road12", "doug's office", "culvert 0", "artesia", "paseo del rio", "leasburg dam", "la8", "n. caballo lake", "n. caballo lake 2", "fish well", "lm bar, 1443 tank - poss south seco tank", "1342 tank - possibly south seco tank", "pague well", "n. seco well", "sissel tank", "polo tank", "n. caballo lake 3", "percha dam state park", "southwell tank", "unid jornada tank", "frog tank"), lat = c(41.569722, 41.691705, 41.593133, 41.665167, 41.672974, 41.621208, 41.650494, 41.638883, 41.665733, 41.691276, 41.643053, 41.712117, 35.35072, 35.343463, 35.40982, 35.42517, 35.41078, 44.507673, 44.515034, 44.498724, 44.627464, 44.496218, 44.493559, 44.400566, 44.376923, 31.700297, 31.484564, 31.46769, 31.444947, 31.320851, 31.255013, NA, NA, 32.739025, 33.090787, 32.355656, NA, 33.499471, 33.495826, 33.493784, 33.133323, 33.337197, 33.495381, 33.331613, 33.164875, 33.194709, NA, 33.169416, 33.507536, NA, 32.991219, 33.050854, 33.151993, 32.487049, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), long = c(-80.4525, -80.50082, -80.356283, -80.51465, -80.513485, -80.469131, -80.423694, -80.495117, -80.511781, -80.512482, -80.430097, -80.389028, -86.13993, -86.129624, -86.06988, -86.07452, -86.07966, -73.162141, -73.162287, -73.171326, -73.237832, -72.890135, -73.242776, -73.23624, -73.161736, -92.970663, -93.174226, -93.147099, -93.093505, -93.134511, -93.24072, NA, NA, -106.71321, -107.566088, -106.404037, NA, -106.359682, -106.293403, -106.327882, -106.505509, -106.313692, -106.292815, -107.64518, -107.61595, -106.696956, NA, -107.684043, -106.205415, NA, -106.519744, -107.529161, -107.203096, -106.924443, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), site_comments = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Radium Springs Exit", "Jornada", "Ladder Ranch", "White Sand Missile Range", "Ladder Ranch", "White Sand Missile Range", "White Sand Missile Range", "White Sand Missile Range", "White Sand Missile Range", "White Sand Missile Range", "White Sand Missile Range", "Ladder Ranch", "Ladder Ranch", "White Sand Missile Range", "White Sand Missile Range", "Ladder Ranch", "White Sand Missile Range", "Leesburg Dam", "White Sand Missile Range", "Ladder Ranch", "Elephant Butte", "Leasburg Dam", NA, "Williamsburg", "Williamsburg", "Ladder Ranch", "Ladder Ranch", "Ladder Ranch", "Ladder Ranch", "Ladder Ranch", "Ladder Ranch", "Ladder Ranch", "Williamsburg", "Caballo Lake", "Jornada", "Jornada", "Jornada"), location = c("usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa", "usa"), utm_zone = c(17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 16, 16, 16, 16, 16, 18, 18, 18, 18, 18, 18, 18, 18, 15, 15, 15, 15, 15, 15, NA, NA, 13, 13, 13, NA, 13, 13, 13, 13, 13, 13, 13, 13, 13, NA, 13, 13, NA, 13, 13, 13, 13, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)), row.names = c(NA, -69L), class = c("tbl_df", "tbl", "data.frame"))

现有单区转换脚本 #

v <- vect(data, c("long", "lat"), crs="+proj=longlat")
u <- terra::project(v, " +proj=utm +zone=23")
utm <- crds(u)

i <- which(!is.na(data$long))
data[i, c("utme", "utmn")] <- utm

data <- data %>% 
  select(!c(long, lat))

批量转换解决方案 #

方法一：基础循环实现

逻辑清晰，适合新手理解：

# 加载依赖包
library(terra)
library(dplyr)

# 初始化UTM坐标列
data$utme <- NA_real_
data$utmn <- NA_real_

# 获取所有有效UTM带（去重并排除NA）
valid_zones <- unique(na.omit(data$utm_zone))

# 遍历每个UTM带进行转换
for (zone in valid_zones) {
  # 筛选当前UTM带且经纬度非空的行索引
  target_rows <- which(data$utm_zone == zone & !is.na(data$lat) & !is.na(data$long))
  
  # 如果没有符合条件的行，跳过当前循环
  if (length(target_rows) == 0) next
  
  # 转换为terra矢量对象
  point_vect <- vect(data[target_rows, ], geom = c("long", "lat"), crs = "+proj=longlat +datum=WGS84")
  
  # 投影到对应UTM带
  utm_vect <- terra::project(point_vect, paste0("+proj=utm +zone=", zone, " +datum=WGS84"))
  
  # 提取UTM坐标并赋值回原数据框
  utm_coords <- crds(utm_vect)
  data[target_rows, c("utme", "utmn")] <- utm_coords
}

# 可选：移除原经纬度列
data <- data %>% select(-long, -lat)

方法二：dplyr分组处理（更简洁高效）

适合熟悉tidyverse语法的用户：

library(terra)
library(dplyr)
library(purrr)

# 处理数据并转换坐标
data_processed <- data %>%
  # 保留有效数据行
  filter(!is.na(utm_zone), !is.na(lat), !is.na(long)) %>%
  # 按UTM带分组处理
  group_by(utm_zone) %>%
  mutate(
    # 每组内完成坐标转换
    utm_data = list({
      v <- vect(cur_data(), geom = c("long", "lat"), crs = "+proj=longlat +datum=WGS84")
      u <- terra::project(v, paste0("+proj=utm +zone=", first(utm_zone), " +datum=WGS84"))
      crds(u)
    }),
    # 拆分坐标列
    utme = map_dbl(utm_data, ~ .[1]),
    utmn = map_dbl(utm_data, ~ .[2])
  ) %>%
  ungroup() %>%
  # 合并回原数据（保留含NA的行）
  right_join(data, by = names(data)) %>%
  # 清理临时列和原经纬度列
  select(-long, -lat, -utm_data)

注意事项 #

• 转换时指定+datum=WGS84确保坐标基准统一，避免误差
• 两种方法都会自动跳过utm_zone、lat或long为NA的行，防止报错
• 若处理南半球数据，需在投影字符串中添加+south参数（例如+proj=utm +zone=35 +south +datum=WGS84）

内容的提问来源于stack exchange，提问作者Eizy