如何在R中实现极地立体投影与经纬度互转及海冰数据提取？

问题概述 #

需从南大洋指定区域（CCAMLR 88.1/88.2/88.3区）提取月度合成海冰浓度数据，目标数据集采用极地立体投影（polar stereographic projection），但项目全程使用经纬度（lat-lon）坐标系，需完成两项核心转换：

• 将经纬度定义的提取多边形转换为极地投影坐标，用于数据子集提取
• 将提取后的极地投影数据的x/y坐标转回经纬度，用于后续绘图

此前仅找到Python实现方案，现需R语言的便捷解决方法。

核心解决方案：使用sf包处理投影转换 #

sf是R语言处理空间数据的标准工具，支持任意投影间的转换。针对南极极地立体投影，标准EPSG编码为3031（WGS 84 / Antarctic Polar Stereographic），经纬度坐标系采用标准EPSG:4326（WGS84）。

步骤1：安装并加载必要包 #

install.packages(c("sf", "ggplot2", "dplyr"))
library(sf)
library(ggplot2)
library(dplyr)

步骤2：经纬度多边形转极地投影 #

将你定义的CCAMLR区域多边形转换为sf对象，再转换到极地投影：

# 定义经纬度多边形（原代码）
Area_88.1 <- data.frame(
  long = c(150, 190, 190, 150),
  lat = c(-60, -60, -85, -85))
Area_88.2 <- data.frame(
  long = c(-170, -105, -105, -170),
  lat = c(-60, -60, -85, -85))
Area_88.3 <- data.frame(
  long = c(-105, -70, -70, -105),
  lat = c(-60, -60, -76, -76))

# 转换为sf对象，设置经纬度CRS（EPSG:4326）
area_88.1_sf <- st_as_sf(Area_88.1, coords = c("long", "lat"), crs = 4326) %>% 
  st_polygonize()
area_88.2_sf <- st_as_sf(Area_88.2, coords = c("long", "lat"), crs = 4326) %>% st_polygonize()
area_88.3_sf <- st_as_sf(Area_88.3, coords = c("long", "lat"), crs = 4326) %>% st_polygonize()

# 转换为南极极地立体投影（EPSG:3031）
area_88.1_polar <- st_transform(area_88.1_sf, crs = 3031)
area_88.2_polar <- st_transform(area_88.2_sf, crs = 3031)
area_88.3_polar <- st_transform(area_88.3_sf, crs = 3031)

转换后的area_88.x_polar可直接用于提取极地投影下的海冰数据子集。

步骤3：提取后的数据转回经纬度 #

假设提取到的海冰数据包含x（极地投影x坐标）、y（极地投影y坐标）、sea_ice_conc（海冰浓度）列，转换为经纬度：

# 示例提取后的数据结构
sea_ice_data <- data.frame(
  x = c(-1000000, -800000),  # 极地投影x
  y = c(-2000000, -1800000), # 极地投影y
  sea_ice_conc = c(0.6, 0.8) # 海冰浓度
)

# 转换为sf点对象，设置极地投影CRS
sea_ice_sf <- st_as_sf(sea_ice_data, coords = c("x", "y"), crs = 3031)

# 转回经纬度坐标系
sea_ice_latlon <- st_transform(sea_ice_sf, crs = 4326)

# 提取lat/lon列到普通数据框
sea_ice_latlon_df <- sea_ice_latlon %>% 
  mutate(
    lon = st_coordinates(.)[,1],
    lat = st_coordinates(.)[,2]
  ) %>% 
  st_drop_geometry()

处理后的sea_ice_latlon_df可直接用于经纬度绘图流程。

步骤4：整合到绘图流程 #

用sf的geom_sf简化投影绘图，替换原ggplot代码中的geom_polygon：

# 加载世界地图的sf对象
wm_sf <- st_as_sf(map_data("world"), coords = c("long", "lat"), crs = 4326) %>% 
  st_set_geometry(st_combine(st_geometry(.))) %>% 
  st_cast("POLYGON")

# 绘图
ggplot() +
  # 世界陆地
  geom_sf(data = wm_sf, fill = "grey", colour = "black", alpha = 1) +
  # CCAMLR区域（经纬度sf对象直接绘制）
  geom_sf(data = area_88.1_sf, fill = "blue", colour = "black", alpha = 0.5) +
  geom_sf(data = area_88.2_sf, fill = "red", colour = "black", alpha = 0.5) +
  geom_sf(data = area_88.3_sf, fill = "green", colour = "black", alpha = 0.5) +
  # 设置正交投影（匹配原绘图效果）
  coord_sf(crs = st_crs("+proj=ortho +lat_0=-90 +lon_0=210"), xlim = c(-60, -240), ylim = c(-90, -56)) +
  # 样式设置
  theme(
    panel.background = element_blank(),
    panel.grid.major = element_line(linewidth = 0.25, linetype = "dashed", colour = "black"),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank()
  )

额外说明 #

若数据集使用非标准极地投影参数，可手动定义投影字符串，例如NSIDC旧版极地投影：

nsIDC_polar_crs <- st_crs("+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

用该CRS替代之前的EPSG:3031即可适配。

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