如何基于R语言对已栅格化的多边形站点按Site_Type分类？

嘿，你已经完成了从CSV到Robinson投影多边形的大部分工作，核心问题就是在转换过程中丢失了原始的站点类型（Site_Type）属性——现在sp_robinson里只有layer字段，没有咱们需要的site信息。下面是分步的修正方案，帮你把类型属性找回来：

一、先创建带完整属性的缓冲区多边形 #

你之前创建nodes的时候，只是单纯做了缓冲区操作，没保留原始数据的属性。咱们先把CSV转成sf对象，再创建缓冲区，这样就能把site、radius这些字段都带过来：

# 把原始数据转成sf空间对象，指定经纬度坐标和CRS
sp_sf <- st_as_sf(sp_csv_data, coords = c("lon", "lat"), crs = longlat)
# 给每个站点创建缓冲区，同时保留所有原始属性
nodes <- st_buffer(sp_sf, dist = sp_sf$radius)

现在nodes里就有site字段了，接下来栅格化的时候就能把这个属性关联上。

二、处理单/多站点类型的栅格化 #

你的数据里有NES;BRE这种多类型的站点，这里提供两种处理方式，按需选择：

方式1：拆分多类型为独立多边形 #

如果需要把一个多类型站点拆成多个分别对应不同类型的多边形（比如把NES;BRE拆成两个多边形，一个标NES，一个标BRE），可以先拆分site字段：

library(tidyr) # 需要用到unnest函数
library(stringr) # 需要用到str_split函数

# 拆分多类型字符串，拆成单独的行
sp_sf_split <- sp_sf %>%
  mutate(site = str_split(site, ";")) %>% # 按分号拆分类型
  unnest(site) # 把拆分后的类型展开成独立行

# 基于拆分后的数据创建缓冲区，每个类型对应一个多边形
nodes_split <- st_buffer(sp_sf_split, dist = sp_sf_split$radius)

方式2：保留原始多类型字符串 #

如果想保留NES;BRE这种原始的多类型格式，咱们可以给每个类型组合分配唯一ID，栅格化后再映射回类型名称：

# 给每个类型组合分配唯一整数ID（方便栅格化处理）
sp_sf$site_id <- as.factor(sp_sf$site) %>% as.integer()
# 创建带ID的缓冲区
nodes <- st_buffer(sp_sf, dist = sp_sf$radius)

# 栅格化时指定用site_id作为字段，这样栅格就会带上类型ID
sp_raster <- rasterize(nodes, rs, field = "site_id")

# 把栅格转成sf多边形，并关联回原始类型名称
site_map <- distinct(sp_sf, site_id, site) # 建立ID和类型的映射表
species_sp <- as(sp_raster, "SpatialPolygonsDataFrame") %>%
  st_as_sf() %>%
  left_join(site_map, by = c("layer" = "site_id")) # 把ID替换成类型名称

三、修正后续投影与日界线处理 #

把上面的代码替换你原来的Creating the nodes到converting into an sf spatial polygon dataframe部分，之后的投影、日界线切分代码可以直接保留。最后得到的sp_robinson就会包含site字段啦！

完整的修正后代码片段：

# 替换原有的Creating the nodes到converting部分
sp_sf <- st_as_sf(sp_csv_data, coords = c("lon", "lat"), crs = longlat)
sp_sf$site_id <- as.factor(sp_sf$site) %>% as.integer()
nodes <- st_buffer(sp_sf, dist = sp_sf$radius)

# 栅格化
rs <- raster(ncol = 360*2, nrow = 180*2)
rs[] <- 1:ncell(rs)
crs(rs) <- CRS(longlat)
sp_raster <- rasterize(nodes, rs, field = "site_id")
sp_raster <- resample(sp_raster, rs, resample = "ngb")

# 转成sf并关联类型
site_map <- distinct(sp_sf, site_id, site)
species_sp <- as(sp_raster, "SpatialPolygonsDataFrame") %>%
  st_as_sf() %>%
  left_join(site_map, by = c("layer" = "site_id")) %>%
  st_set_crs(longlat)

# 以下是你原来的投影和日界线处理代码，直接保留
polygon <- st_polygon(x = list(rbind(c(-0.0001, 90), c(0, 90), c(0, -90), c(-0.0001, -90), c(-0.0001, 90)))) %>% st_sfc() %>% st_set_crs(longlat)
sp_robinson <- species_sp %>% st_difference(polygon) %>% st_transform(crs = rob_pacific)

# 修复南极洲分割线
bbox1 <- st_bbox(sp_robinson)
bbox1[c(1,3)] <- c(-1e-5,1e-5)
polygon1 <- st_as_sfc(bbox1)
crosses1 <- sp_robinson %>% st_intersects(polygon1) %>% sapply(length) %>% as.logical %>% which
sp_robinson[crosses1, ] %<>% st_buffer(0)

四、验证结果 #

现在你可以用head(sp_robinson)查看属性，会发现site字段已经存在；也可以用plot(sp_robinson["site"])直接按站点类型可视化，每个类型会用不同颜色区分。

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