鲨类家域MCP坐标转换异常：UTM转WGS84适配谷歌地图故障排查

问题描述 #

我正在用最小凸多边形（MCP）研究鲨类家域，已成功生成MCP，但叠加谷歌地图时坐标完全错误。问题出在**UTM 40区转WGS84（EPSG:4326）**的步骤：原本应位于查戈斯群岛（坐标约71.96, -5.24），转换后变成52.51, -4.72e-05，完全偏离目标区域。

相关代码 #

转换为UTM 40N（EPSG:32640） #

此步骤输出坐标正常：

SF32map_cat2 <- st_as_sf(SP_cat2, crs = 32640) 
SF32map_cat2MCP<- st_as_sf(SP_cat2MCP, crs = 32640)

转换为WGS84（EPSG:4326） #

执行此步骤后坐标异常：

SFmap_cat2 <- st_transform(SF32map_cat2, crs = 4326) 
SFmap_cat2MCP <- st_transform(SF32map_cat2MCP, crs = 4326)

绘图代码 #

ggmap(chagos_map) + 
  geom_sf(data = SFmap_cat2, aes(color = as.factor(code)), size = 2, inherit.aes = FALSE) +
  geom_sf(data = SFmap_cat2MCP, fill = NA, color = "black", alpha = 0.5, inherit.aes = FALSE) +
  ggtitle("MCP Silvertip Cat 2 Chagos") +
  theme_minimal()

已排查情况 #

• 单独绘制MCP形状正确，但坐标错误
• 确认数据类型与结构无问题
• 尝试改用st_as_sf()直接转换，结果仍坐标错误

数据示例 #

SP_cat2（SpatialPointsDataFrame） #

new("SpatialPointsDataFrame", data = structure(list(code =     c("12964", 
"12964", "12964", "12964", "12964", "12964")), row.names =     c(NA, 
6L), class = "data.frame"), coords.nrs = 2:3, coords =     structure(c(71.9684, 
71.9684, 71.9684, 71.9684, 71.9684, 71.9698, -5.2417,     -5.2417, 
-5.2417, -5.2417, -5.2417, -5.3302), dim = c(6L, 2L),     dimnames = list(
NULL, c("receiver_lon", "receiver_lat"))), bbox =     structure(c(71.9684, 
-5.3302, 71.9698, -5.2417), dim = c(2L, 2L), dimnames = list(
c("receiver_lon", "receiver_lat"), c("min", "max"))),     proj4string = new("CRS", 
projargs = "+proj=utm +zone=40 +datum=WGS84 +units=m         +no_defs"))

SP32map_cat2MCP（sf对象） #

structure(list(code = c("12964", "12964", "12964", "12964", "12964", 
"12964"), geometry = structure(list(structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9698, -5.3302
), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 71.9684, ymin = -5.3302, 
xmax = 71.9698, ymax = -5.2417), class = "bbox"), crs = structure(list(
input = "+proj=utm +zone=40 +datum=WGS84 +units=m +no_defs", 
wkt = "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6326]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"UTM zone 40N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",57,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]],\n        ID[\"EPSG\",16040]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(code = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = c(NA, 
6L), class = c("sf", "data.frame"))

解决方案 #

问题根源是你把地理坐标（经纬度）当成UTM坐标赋值给了EPSG:32640，且选错了UTM带的EPSG代码。

核心错误点 #

1. 查戈斯群岛位于南半球，对应的UTM带是40S，EPSG代码为32740（32640是北半球的40N）。
2. 原始数据的receiver_lon和receiver_lat是经纬度（WGS84），单位是度，不是UTM的米制坐标，不能直接以UTM坐标系导入。

正确操作步骤 #

1. 先以WGS84（EPSG:4326）导入原始数据：

# 把SP_cat2转为WGS84坐标系的sf对象
SF_wgs84_cat2 <- st_as_sf(SP_cat2, crs = 4326)

2. 转换到正确的UTM带（32740）生成MCP：

library(adehabitatHR)
# 转UTM 40S
SF_utm_cat2 <- st_transform(SF_wgs84_cat2, crs = 32740)
# 生成MCP
SP_cat2MCP <- mcp(SF_utm_cat2, percent = 95)
# 转回WGS84用于绘图
SF_mcp_wgs84 <- st_as_sf(SP_cat2MCP) %>% st_transform(crs = 4326)

3. 重新获取底图并绘图：

# 确保底图是WGS84坐标
chagos_map <- get_map(location = c(lon = 71.96, lat = -5.24), zoom = 10, maptype = "terrain")
# 叠加绘制
ggmap(chagos_map) + 
  geom_sf(data = SF_wgs84_cat2, aes(color = as.factor(code)), size = 2, inherit.aes = FALSE) +
  geom_sf(data = SF_mcp_wgs84, fill = NA, color = "black", alpha = 0.5, inherit.aes = FALSE) +
  ggtitle("MCP Silvertip Cat 2 Chagos") +
  theme_minimal()

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