Combine units in a shapefile while keeping others granular

I have a 5-digit postcode shapefile for Germany. Big number 1-digits postcodes are similar to German states. I read shapefile data with rgdal thus having a SpatialPolygonsDataFrame. I have only data for some part of Germany, i.e. some postcodes. The data I have I like to show on a granular 5-digit level. Using leaflet to create a map, it makes very long time for me, to plot all of the almost 10.000 postcodes. Thus I like to "summarize"/"combine"/"merge" the outer border of those postcodes where I don't have data (where the value is NA).

# German postcode shapes 

# Create temp files
temp <- tempfile()
temp2 <- tempfile()

# Download the zip file and save to 'temp' 
URL <- "https://downloads.suche-postleitzahl.org/v2/public/plz-5stellig.shp.zip"
download.file(URL, temp)

# Unzip the contents of the temp and save unzipped content in 'temp2'
unzip(zipfile = temp, exdir = temp2)

# Read shape file 
library(rgdal)
GER_postcode <- readOGR(temp2)

head(GER_postcode@data$note)

# Create subsample 
library(tidyverse)

GER_postcode@data$einwohner2 <- ifelse(substr(GER_postcode@data$plz, 1, 1) %in% c("0", "1", "7"), GER_postcode@data$einwohner, NA)

# Plot Subsample 
library(leaflet)

qpal <- colorBin("Reds", GER_postcode@data$einwohner2, bins=10)

leaflet(GER_postcode) %>%
  addPolygons(stroke = TRUE,opacity = 1,fillOpacity = 0.5, smoothFactor = 0.5,
              color="black",fillColor = ~qpal(einwohner2),weight = 1) %>%
  addLegend(values=~einwohner2,pal=qpal,title="Population")

How can I make the map show those postcode shapes with values and combine all other where the value is NA?

I was looking at library(rgeos) and gUnaryUnion() that units all units in a shapefile to outer borders. Although I only need to do this on a subset.

I personally would use filter / subset to drop irrelevant polygons and use country outline as a separate background layer. Outline can be built by making an union of postcode shapes, though it's likely faster to download one (or install some package that bundles country shapes)

As rgdal and regos will be retired in just few months, the example below uses sf (and optional geos package ) instead.

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(dplyr)
library(leaflet)

URL <- "https://downloads.suche-postleitzahl.org/v2/public/plz-5stellig.shp.zip"

# use GDAL virtual file systems to load zipped shapefile from remote url
GER_postcode <- paste0("/vsizip//vsicurl/", URL) %>%  read_sf()

# country outline from giscoR
GER_outline <- giscoR::gisco_get_countries(country = "DE")

# if it must be built from input shapes, we can use sf::st_union() or switch  
# to geos (the R library), it performs bit better in this case:

# GER_outline <- as_geos_geometry(GER_postcode) %>% 
#   geos_make_collection() %>% 
#   geos_unary_union() %>% 
#   st_as_sfc()

GER_postcode_subsample <- GER_postcode %>% filter(substr(plz, 1, 1) %in% c(0, 1, 7))

# Plot Subsample 
qpal <- colorBin("Reds", GER_postcode_subsample$einwohner, bins=10)

leaflet(GER_postcode_subsample) %>%
  addPolygons(data = GER_outline, smoothFactor = 0.5,
              color="black", fillColor = "#808080", weight = 1) %>% 
  addPolygons(stroke = TRUE,opacity = 1,fillOpacity = 0.5, smoothFactor = 0.5,
              color="black",fillColor = ~qpal(einwohner),weight = .2) %>%
  # append NA to values to include it in legend
  addLegend(values=c(~einwohner, NA),pal=qpal,title="Population")

^{Created on 2023-07-17 with reprex v2.0.2}