R提取南极海冰边缘轮廓转多边形并导入QGIS的问题求助

问题背景 #

下载了南极海冰浓度GeoTiff文件，需提取15%浓度对应的海冰边缘轮廓，转换为可在QGIS中打开并重投影的格式。原代码尝试将轮廓线转为SpatialPointsDataFrame后直接转SpatialPolygons，触发错误：

Error in as(ie.pt, "SpatialPolygons"): no method or default for coercing “SpatialPointsDataFrame” to “SpatialPolygons”

错误原因 #

SpatialPointsDataFrame是无序点集，丢失了原轮廓线的拓扑连接信息，无法直接转换为多边形。必须基于原始轮廓线的分组结构，将同一条线的点按顺序组成闭合环后再构建多边形。

解决方案 #

方法1：基于sp包修正原工作流 #

保留原始轮廓线的分组信息，将每条线的点按顺序转换为多边形：

library(raster)
library(sp)
library(sf)

# 加载并预处理海冰浓度数据
sic <- raster('Environmental_Data/SIC/AMSR2/asi-AMSR2-s3125-20220107-v5.4.tif')
sic[sic > 100] <- NA  # 标记陆地为NA
sic <- crop(sic, extent(sic)*c(0.5,0.5,0,1))  # 裁剪到目标区域

# 提取15%浓度的海冰边缘轮廓线
ie <- rasterToContour(sic, levels=15)

# 将每条轮廓线转换为多边形
polygons_list <- lapply(1:length(ie), function(i) {
  # 获取当前线的坐标点
  line_coords <- ie@lines[[i]]@Lines[[1]]@coords
  # 构建Polygon（确保闭合，若线已闭合可省略rbind步骤）
  poly <- Polygon(rbind(line_coords, line_coords[1,]))
  Polygons(list(poly), ID=as.character(i))
})

# 转为SpatialPolygons并添加属性
ie_poly <- SpatialPolygons(polygons_list)
ie_poly_df <- SpatialPolygonsDataFrame(ie_poly, data=data.frame(id=1:length(ie), row.names=1:length(ie)))

# 导出为Shapefile（QGIS可直接打开）
writeOGR(ie_poly_df, dsn="sea_ice_edge_polygons", layer="ice_edge", driver="ESRI Shapefile", overwrite_layer=TRUE)

方法2：使用sf包简化流程（推荐） #

sf包处理空间数据更直观，支持直接将轮廓线转为多边形并导出：

library(raster)
library(sf)

# 加载并预处理数据（同前）
sic <- raster('Environmental_Data/SIC/AMSR2/asi-AMSR2-s3125-20220107-v5.4.tif')
sic[sic > 100] <- NA
sic <- crop(sic, extent(sic)*c(0.5,0.5,0,1))

# 提取轮廓线并转为sf对象
ie_sf <- st_as_sf(rasterToContour(sic, levels=15))

# 将线转为多边形（需确保线是闭合环，若有不闭合情况可先处理）
ie_poly_sf <- st_polygonize(ie_sf)

# 导出为GeoPackage（QGIS更推荐的格式），或Shapefile
st_write(ie_poly_sf, "sea_ice_edge_polygons.gpkg", layer="ice_edge", delete_layer=TRUE)
# 若需Shapefile：
# st_write(ie_poly_sf, "sea_ice_edge_polygons/ice_edge.shp", delete_layer=TRUE)

补充说明 #

• 导出的文件可直接拖入QGIS，在QGIS中通过「投影变换」工具完成重投影操作
• 若轮廓线存在不闭合情况，可在sf中使用st_close函数先闭合线条再转多边形

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