r - 在 R 中加入不同的 SpatialPolygonsDataFrame 对象时出现问题

Question

我有一个西类牙北部城镇的形状文件，我必须加入组(西类牙语中的自治市或 comarcas)。我已经使用 sf 包中的 st_union 成功加入了它们(每个都是它们自己的具有单个多边形的 SpatialPolygonsDataFrame 对象)。我分别绘制了每个城市，它们看起来不错。

但是，一旦我想将城市合并到一个具有多个多边形的 SpatialPolygonsDataFrame 对象中，我就没法做到了。我已经尝试了三种主要基于这个答案的方法:https://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r还有这个https://gis.stackexchange.com/questions/141469/how-to-convert-a-spatialpolygon-to-a-spatialpolygonsdataframe-and-add-a-column-t– 如果我使用 raster::union 它会抛出错误.rowNamesDF<- 中的错误(x, value = value) : 无效的 'row.names' 长度– 如果我使用简单的 rbind，它会抛出错误SpatialPolygonsDataFrame(pl, df, match.ID = FALSE) 错误:对象长度不匹配:pl 有 7 个多边形对象，但 df 有 4 行或 6/11 个城市的类似情况。– 如果我尝试一种 lapply 方法(更复杂)，它似乎有效，但我使用传单绘制了它，在尝试 raster::union 或 rbind 时出现错误的市政当局看起来不应该/不看起来当我单独绘制它们时，它们会这样做。

** 市政府 1 和 2 工作正常。例如 3 和 4 没有。 **

下面是重现我的代码所需的两个文件的链接:– 形状文件链接:https://www.dropbox.com/sh/z9632hworbbchn5/AAAiyq3f_52azB4oFeU46D5Qa?dl=0– 链接到包含从城镇到自治市的映射的 xls 文件:https://www.dropbox.com/s/4w3fx6neo4t1l3d/listado-comarcas-gipuzkoa.xls?dl=0

还有我的代码:

library(tidyverse)
library(magrittr)
library(sf)
library(ggplot2)
library(lwgeom)
library(readxl)
library(raster)

#Read shapefile
mapa_municip <- readOGR(dsn = "UDALERRIAK_MUNICIPIOS/UDALERRIAK_MUNICIPIOS.shp")
mapa_municip <- spTransform(mapa_municip, CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'))
mapa_municip <- st_as_sf(mapa_municip)

#Read excel that contains mapping from town to municioalities
muni2com     <- read_excel("listado-comarcas-gipuzkoa.xls",
                           sheet=1,
                           range="A1:C91", 
                           col_names = T)

comarcas <- list()

count <- 0

for (i in unique(muni2com$Comarca)[1:4]){
  
  count <- count + 1 
  
  for (k in unique(muni2com$Municipios[muni2com$Comarca==i])){
    
    if (k == unique(muni2com$Municipios[muni2com$Comarca==i])[1]){ # if 1st case, keep this town
      temp <- mapa_municip[mapa_municip$MUNICIPIO==k,]
    } 
    if (k != unique(muni2com$Municipios[muni2com$Comarca==i])[1]){ # otherwise, join w previous ones
      temp <- sf::st_union(temp, mapa_municip[mapa_municip$MUNICIPIO==k,])
      
    }
    
  }
  
  comarcas[[count]] <-  spTransform(as(temp, "Spatial"), CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'))
  comarcas[[count]]@data <- data.frame(comarca = i)
    
}

IDs <- sapply(comarcas, function(x)
  slot(slot(x, "polygons")[[1]], "ID"))

#Checking
length(unique(IDs)) == length(comarcas)

dfIDs <- data.frame(comarca = IDs)
  
#Making SpatialPolygons from list of polygons
comarcas2 <- SpatialPolygons(lapply(comarcas,
                               function(x) slot(x, "polygons")[[1]]))

# Try to coerce to SpatialPolygonsDataFrame (will throw error)
p.df <- data.frame( comarca = unique(muni2com$Comarca)[1:4])
p <- SpatialPolygonsDataFrame(comarcas2, p.df) 

# Extract polygon ID's
( pid <- sapply(slot(comarcas2, "polygons"), function(x) slot(x, "ID")) )

# Create dataframe with correct rownames
( p.df <- data.frame( comarca = unique(muni2com$Comarca)[1:4], row.names = pid) )    

# Try coertion again and check class
comarcas3 <- SpatialPolygonsDataFrame(comarcas2, p.df)
class(comarcas3) 

#Leaflet map 
leaflet( options = leafletOptions(zoomControl = F,
                                  zoomSnap  = 0.1 ,
                                  zoomDelta = 1
),
data = comarcas3,
) %>%
  addProviderTiles(provider="CartoDB.Positron") %>%
  htmlwidgets::onRender("function(el, x) {
                                       L.control.zoom({ position: 'topright' }).addTo(this)
                                   }") %>%
  clearShapes() %>%
  addPolygons(fillColor = "gray",
              opacity   = 0.8,
              weight = 0.3,
              color = "white",
              fillOpacity = 0.95, 
              smoothFactor = 0.5,
              label = ~comarca,
              highlight = highlightOptions(
                weight = 1.5,
                color = "#333333",
                bringToFront = T),
              layerId = ~comarca
  )

** 请注意，如果您在上面绘制 comarcas[[3]] 或 comarcas[[4]] 而不是 comarcas3，那么这些城市的形状将完全不同。**

我真的很感激你能给我任何提示，我已经解决了好几天了，我无法解决它。我假设问题是由于 rbind 给出的错误引起的，这似乎是最有用的错误，但我不知道它是什么意思。非常感谢您。

Answer 1

您是否绝对需要使用旧的 {sp} 包工作流程？

如果不是，则使用基于纯 {sf} 的工作流将市政当局分解为 comarcas 可能更容易 - 按 comarca 列分组，然后汇总就可以了。

考虑这段代码:

library(tidyverse)
library(sf)
library(readxl)
library(leaflet)

#Read shapefile
mapa_municip <- st_read("UDALERRIAK_MUNICIPIOS.shp") %>% 
  st_transform(4326)


#Read excel that contains mapping from town to municioalities
muni2com <- read_excel("listado-comarcas-gipuzkoa.xls",
                           sheet=1,
                           range="A1:C91", 
                           col_names = T)

# dissolving comarcas using sf / dplyr based workflow
comarcas <- mapa_municip %>% 
  inner_join(muni2com, by = c("MUNICIPIO" = "Municipios")) %>% 
  group_by(Comarca) %>% 
  summarise() %>% # magic! :)))
  ungroup()


leaflet(comarcas) %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(color = "red",
              label = ~ Comarca)