Problem joining different SpatialPolygonsDataFrame objects in R

Question

I have a shape file of towns in the north of Spain that I have to join into groups (municipalities or comarcas in Spanish). I've used st_union from the sf package to join them successfully (and each one is their own SpatialPolygonsDataFrame object with a single polygon). I plot each of the municipalities individually and they look fine.

However, once I want to combine the municipalities into a single SpatialPolygonsDataFrame object with multiple polygons, I can't for the life of me manage to do it. I've tried three approaches mostly based on this answer: https://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r and this one https://gis.stackexchange.com/questions/141469/how-to-convert-a-spatialpolygon-to-a-spatialpolygonsdataframe-and-add-a-column-t – If I use raster::union it throws out the error Error in .rowNamesDF<-(x, value = value) : invalid 'row.names' length – If I use a simple rbind it throws out the error Error in SpatialPolygonsDataFrame(pl, df, match.ID = FALSE) : Object length mismatch: pl has 7 Polygons objects, but df has 4 rows Or something similar for 6/11 of the municipalities. – If I try a lapply approach (more convoluted) it seems to work but one I plot it using leaflet the municipalities that gave the error when trying to raster::union or rbind don't look as they should/don't look as they do when I plot them individually.

** Municipalities 1 and 2 work fine. 3 and 4 for example do not. **

Here's a link to the two files needed to reproduce my code below: – Link to shape files: https://www.dropbox.com/sh/z9632hworbbchn5/AAAiyq3f_52azB4oFeU46D5Qa?dl=0 – Link to xls file that contains the mapping from towns to municipalities: https://www.dropbox.com/s/4w3fx6neo4t1l3d/listado-comarcas-gipuzkoa.xls?dl=0

And my code:

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
  )

** Note how if you plot comarcas[[3]] or comarcas[[4]] above instead of comarcas3 the shape of those municipalities is completely different.**

I'd really appreciate any tips you can give me, I've been at it for days and I can't solve it. I assume the problem is due to the error given by the rbind, which seems to be the most informative one, but I don't know what it means. Thank you very much in advance.

Answer 1

Are you absolutely positively required to use the older {sp} package workflow?

If not it may be easier to dissolve the municipalities into comarcas using a pure {sf} based workflow - grouping by a comarca column, and then summarising will do the trick.

Consider this code:

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)