4

I use to plot maps of Belgium available on GADM (select Belgium) using R.

I import and plot the outside boundary of Belgium using :

belgium <-readRDS("gadm36_BEL_0_sp.rds")
plot(belgium)

Which gives me :

Map of Belgium

I import and plot the boundaries of provinces of Belgium using :

provinces <-readRDS("gadm36_BEL_2_sp.rds")
plot(provinces)

Which gives me :

Map of provinces of Belgium

What i'm trying to have is a dataframe with boundaries of provinces that are NOT outside boundaries of Belgium :

enter image description here

I tried using over(), intersect(), etc but did not founded yet a method to do that. Several approach can be used I guess :

  • substract belgium boundary to the province dataset;
  • Work only in the provinces dataset and keep only commons boundaries;
  • else?

Thanks if you have a solution. Grégoire

massisenergy
  • 1,764
  • 3
  • 14
  • 25
Gregoire Vincke
  • 313
  • 2
  • 9

1 Answers1

3

I downloaded the sf formatted files from that site (https://www.gadm.org/download_country_v3.html), since the sf package is a bit easier to deal with.

library(dplyr)
library(sf)

provinces <- readRDS("gadm36_BEL_2_sf.rds")

interiors <- st_intersection(provinces) %>% 
  filter(n.overlaps > 1)

interiors

# Number of columns truncated for clarity:
#   interiors %>% select(VARNAME_2, geometry, n.overlaps)

Simple feature collection with 30 features and 2 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 2.851679 ymin: 49.8004 xmax: 6.033082 ymax: 51.35568
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
                                                                                                            VARNAME_2
1                                                                            Amberes|Antuérpia|Antwerp|Anvers|Anversa
2                                                                            Amberes|Antuérpia|Antwerp|Anvers|Anversa
3  Brussel Hoofstadt|Brusselse Hoofdstedelijke Gewest|Brüssel|Bruxelas|Région de Bruxelles-Capitale|Brussels|Bruselas
4                                                                                                   Limbourg|Limburgo
5                   Flandres Oriental|Fiandra Orientale|Flandes Oriental|Flandre orientale|East Flanders|Ost Flandern
6                                                                            Amberes|Antuérpia|Antwerp|Anvers|Anversa
7                                                                            Amberes|Antuérpia|Antwerp|Anvers|Anversa
8                                                                            Amberes|Antuérpia|Antwerp|Anvers|Anversa
9                   Flandres Oriental|Fiandra Orientale|Flandes Oriental|Flandre orientale|East Flanders|Ost Flandern
10                                                Brabant Flamand|Brabante Flamenco|Brabante Flamengo|Flemish Brabant
   n.overlaps                       geometry
1           2 MULTILINESTRING ((5.239571 ...
2           2 MULTILINESTRING ((4.327078 ...
3           2 MULTILINESTRING ((4.403365 ...
4           2 MULTILINESTRING ((5.117446 ...
5           2 MULTILINESTRING ((4.243931 ...
6           3       POINT (4.994605 51.0414)
7           3      POINT (4.243931 51.04332)
8           2 MULTILINESTRING ((4.994605 ...
9           2 MULTILINESTRING ((3.466959 ...
10          2 MULTILINESTRING ((5.025736 ...

To check with a plot:

plot(interiors$geometry)

enter image description here

What you're doing here is looking for the spatial intersection of the provinces with every other province. Then you filter out the intersections where it's just a province overlapping itself (n.overlaps == 1). That way you only get the interior borders where one or more provinces touches another (n.overlaps > 1), but not any province alone (which would be an external border).

This is an updated version of this excellent answer: https://stackoverflow.com/a/47761959/3330437

To remove the circled points (intersections of 3 provinces) in the map and dataset, you can use:

interiors %>% filter(!st_is(., "POINT"))
Brian
  • 7,900
  • 1
  • 27
  • 41
  • 1
    That's just ... perfect ! Thanks ! – Gregoire Vincke Jul 22 '19 at 17:46
  • 1
    Just to complete, for levels 3 (Districts) and 4 (Municipalities), dowloaded at the same adress, i also had to remove MULTIPOINT in addition to POINT. And then, to not having to load sf in addition to sp, i used as(object,"Spatial") to transform the sf object in a sp compliant object, before saving de data in RDS file. – Gregoire Vincke Jul 22 '19 at 18:15