Boundaries

swissBOUNDARIES3D data available from here.

download.file(url = "https://data.geo.admin.ch/ch.swisstopo.swissboundaries3d/swissboundaries3d_2022-01/swissboundaries3d_2022-01_2056_5728.shp.zip",
              destfile = "data-raw/swisstopo/swissboundaries3d_2022-01_2056_5728.shp.zip",
              method = "curl")

unzip("data-raw/swisstopo/swissboundaries3d_2022-01_2056_5728.shp.zip", 
      exdir = "data-raw/swisstopo/swissboundaries3d_2022-01_2056_5728")

unlink("data-raw/swisstopo/swissboundaries3d_2021-01_2056_5728.shp.zip")

# note the source (not raw) to handle #28
st_gg <- 
  st_read("data/swisstopo/swissboundaries3d_2022-01_2056_5728/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_HOHEITSGEBIET.shp",
          as_tibble = TRUE) %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  st_transform(2056) 

Reading layer `swissBOUNDARIES3D_1_3_TLM_HOHEITSGEBIET' from data source 
  `C:\projects\ISPM_excess-mortality-voting\data\swisstopo\swissboundaries3d_2022-01_2056_5728\SHAPEFILE_LV95_LN02\swissBOUNDARIES3D_1_3_TLM_HOHEITSGEBIET.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 2286 features and 23 fields
Geometry type: POLYGON
Dimension:     XYZ
Bounding box:  xmin: 2485410 ymin: 1075268 xmax: 2833858 ymax: 1295934
z_range:       zmin: 193.51 zmax: 4613.686
Projected CRS: CH1903+ / LV95 + LN02 height

st_kt <- 
  st_read("data-raw/swisstopo/swissboundaries3d_2022-01_2056_5728/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp",
          as_tibble = TRUE) %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  st_transform(2056)

Reading layer `swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET' from data source 
  `C:\projects\ISPM_excess-mortality-voting\data-raw\swisstopo\swissboundaries3d_2022-01_2056_5728\SHAPEFILE_LV95_LN02\swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 50 features and 20 fields
Geometry type: POLYGON
Dimension:     XYZ
Bounding box:  xmin: 2485410 ymin: 1075268 xmax: 2833858 ymax: 1295934
z_range:       zmin: 193.51 zmax: 4613.686
Projected CRS: CH1903+ / LV95 + LN02 height

Lakes

Additional dataset of lakes was derived from swissTLM3D_TLM_STEHENDES_GEWAESSER.shp layer from swissTLM3D dataset.

# convert lines to polygons, calculate area

# testing is good here:
# "NAME" = 'Bodensee | Lac de Constance | Lago di Costanza | Lai da Constanza'

st_se <- 
  st_read("data-raw/swisstopo/swisstlm3d_2022-03_2056_5728/TLM_GEWAESSER/swissTLM3D_TLM_STEHENDES_GEWAESSER.shp",
          as_tibble = TRUE) %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  filter(OBJEKTART == "See") %>% 
  group_by(GEW_LAUF_U) %>% 
  summarise(geometry = st_combine(geometry),
            NAME = first(NAME),
            OBJEKTART = first(OBJEKTART)) %>%
  st_cast("MULTILINESTRING") %>% 
  st_polygonize() %>%
  mutate(area = st_area(.)) %>% 
  ungroup() %>% 
  relocate(geometry, .after = last_col())

st_see22_insel <- 
  st_read("data-raw/swisstopo/swisstlm3d_2022-03_2056_5728/TLM_GEWAESSER/swissTLM3D_TLM_STEHENDES_GEWAESSER.shp",
          as_tibble = TRUE) %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  filter(OBJEKTART == "Seeinsel") %>% 
  group_by(GEW_LAUF_U) %>% 
  summarise(geometry = st_combine(geometry),
            NAME = first(NAME),
            OBJEKTART = first(OBJEKTART)) %>%
  st_cast("MULTILINESTRING") %>% 
  st_polygonize() %>% 
  mutate(area = st_area(.)) %>% 
  ungroup() %>% 
  relocate(geometry, .after = last_col())

st_se <- 
  st_difference(st_se, st_union(st_see22_insel)) %>% 
  st_collection_extract("POLYGON") %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  st_transform(2056)

unique(st_geometry_type(st_se))

frq(st_is_valid(st_se))

summary(st_se$area)

plot(st_se[, 2])

rm(st_see22_insel); gc()

Combined

Lakes dataset (selecting only 322 largest features from ThemaKart dataset overlap!) was then used to erase area of municipalities that extends over lakes. Issue is described in #48.

st_gg <- 
  st_difference(st_gg, st_union(st_se_sel)) %>% 
  # small slivers of lakes left
  filter(OBJEKTART != "Kantonsgebiet") %>% 
  st_collection_extract("POLYGON") %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  st_transform(2056) # %>% 
# group_by(BFS_NUMMER) %<>% 
# summarise(BFS_NUMMER = first(BFS_NUMMER),
#           NAME = first(NAME),
#           GEM_TEIL = max(GEM_TEIL)) %>%
# st_cast("MULTIPOLYGON") 

Resolution

swisstopo files (in blue below) are larger in size (and slower to plot) but offer better resolution and will be used for certain spatial operations (like linking to SEP) whereas BfS data (in red) are smaller in size and faster to plot. Example of municipality Bern:

Multipart features

BfS originally stores data as MULTIPOLYGON. In such case municipality split over several disconnected areas is stored as single record. Example of municipality Monthey which consists of four separate areas:

The same municipality in swisstopo data is stored as set of MULTIPOLYGON types and is represented by four separate records instead of one.

That might have some implications in terms of spatial operations (maybe contiguity too?). Data can easily be converted between types if needed.

Disconnected boundaries - Lake Lugano

Two communities in TI Bissone & Melide are disconnected by Lake Lugano in BfS municipality data:

No connectivity

In reality communities are indeed connected:

IRL

That has been manually corrected in GIS and now looks like:

Connected

Importance of this issue is covered here.

Disconnected boundaries - Lake Zurich

Freienbach & Rapperswil-Jona are disconnected on the map of boundaries:

No connectivity

… but connected IRL:

Connectivity

That has been fixed manually in GIS.

Spurious connection

Community of Twann-Tüscherz (in red) is connected to Erlach via small part of the St. Peter’s Island.

St. Peter’s Island

This connection has been removed.

Non-residential municipalities

Excluding municipalities without residents:

• Staatswald Galm
  
• Comunanza Cadenazzo/Monteceneri
  
• Comunanza Capriasca/Lugano

We should not expect any cases there!

ID variables for INLA

gg %<>% 
  mutate(id_space = as.integer(as.factor(GMDNR))) %>% 
  relocate(id_space)

Data sources

More info here.

raum <- read_xlsx("data-raw/BfS/Raumgliederungen.xlsx", 
                  skip = 1) %>% 
  remove_empty(c("rows", "cols")) %>% clean_names() %>% 
  filter(! is.na(bfs_gde_nummer)) %>% 
  select(-bezirksname) %>% 
  rename(GMDNR = bfs_gde_nummer,
         GMDNAME = gemeindename,
         KTNR = kantons_nummer,
         KTNAME = kanton,
         ARGRNR = arbeitsmarktgrossregionen_2018,
         ARNR = arbeitsmarktregionen_2018) %>% 
  select(GMDNR, GMDNAME,
         KTNR, KTNAME, 
         ARGRNR, ARNR) %>% 
  mutate(id_kt = as.integer(as.factor(KTNR))) %>% 
  relocate(id_kt, .after = KTNR) %>% 
  left_join(
    read_xlsx("data-raw/BfS/Raumgliederungen.xlsx", 
              skip = 1, sheet = "CH1+CL_GBAE+2018.0") %>% 
      rename(
        ARGRNR = Code,
        ARGRNAME = Label)
  ) %>% 
  left_join(
    read_xlsx("data-raw/BfS/Raumgliederungen.xlsx", 
              skip = 1, sheet = "CH1+CL_BAE+2018.0") %>% 
      rename(
        ARNR = Code,
        ARNAME = Label) %>% 
      # weird hyphen in some neames?
      mutate(ARNAME = str_replace(ARNAME, fixed("–"), "-"))
  ) %>% 
  relocate(ARGRNAME, .after = ARGRNR) %>% 
  relocate(ARNAME, .after = ARNR) %>% 
  mutate(border = if_else(ARNAME %in% c("Schaffhausen", 
                                        "Delémont", 
                                        "Porrentruy", 
                                        "Saignelégier-Le Noirmont", 
                                        "La Chaux-de-Fonds",
                                        "Le Chenit", 
                                        "Locarno", 
                                        "Lugano") | 
                            ARGRNAME %in% c("Region Basel", "Region Genf"), 
                          1, 0)) 

Border regions

Defined on the basis of this map. Linked by name of the region since some codes seemed to be out of date? (ie. Basel).

Border regions

Number of municipalities that are located within border labour market areas:

border <numeric> 
# total N=2145 valid N=2145 mean=0.20 sd=0.40

Value |    N | Raw % | Valid % | Cum. %
---------------------------------------
    0 | 1725 | 80.42 |   80.42 |  80.42
    1 |  420 | 19.58 |   19.58 | 100.00
 <NA> |    0 |  0.00 |    <NA> |   <NA>

Typology applied to municipalities:

List of municipalities:

gg %<>% left_join(raum %>% select(GMDNR, KTNR, id_kt, KTNAME, 
                                  ARGRNR, ARGRNAME, ARNR, ARNAME,
                                  border)) %>% 
  relocate(geometry, .after = last_col()) 

write_rds(gg, "data/BfS/gg.Rds")

st_write(gg, "data/geo/BfS.gpkg", layer = "st_gg", append = FALSE)

Deleting layer `st_gg' failed
Writing layer `st_gg' to data source `data/geo/BfS.gpkg' using driver `GPKG'
Writing 2145 features with 11 fields and geometry type Multi Polygon.

st_gg %<>% left_join(raum %>% select(GMDNR, KTNAME, 
                                     ARGRNR, ARGRNAME, ARNR, ARNAME,
                                     border)) %>% 
  relocate(geometry, .after = last_col())

write_rds(st_gg, "data/swisstopo/st_gg.Rds")

st_write(st_gg, "data/geo/swisstopo.gpkg", layer = "st_gg", append = FALSE)

Deleting layer `st_gg' failed
Writing layer `st_gg' to data source `data/geo/swisstopo.gpkg' using driver `GPKG'
Writing 2239 features with 10 fields and geometry type Multi Polygon.

Weights

gg_wm_q <- poly2nb(gg, queen = TRUE, row.names = gg$id_space)

gg_wm_q

Neighbour list object:
Number of regions: 2145 
Number of nonzero links: 12186 
Percentage nonzero weights: 0.264854 
Average number of links: 5.681119 

summary(gg_wm_q)

Neighbour list object:
Number of regions: 2145 
Number of nonzero links: 12186 
Percentage nonzero weights: 0.264854 
Average number of links: 5.681119 
Link number distribution:

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
 14  43 187 374 498 409 274 162  95  45  18   7   4   6   2   4   1   1   1 
14 least connected regions:
757 759 820 896 911 992 1002 1118 1155 1408 1558 1582 2072 2127 with 1 link
1 most connected region:
1704 with 19 links

write_rds(gg_wm_q, "data/nb/gg_wm_q.Rds")
nb2INLA("data/nb/gg_wm_q.adj", gg_wm_q)

Region with one neighbour

2271 Meyriez

Region with nineteen neighbours (max case)

5586 Lausanne

Multipart features

Example of 418 Oberburg. Location:

Thanks to using multipart features, connectivity is captured correctly:

Connectivity over Zurichsee

1322 Freienbach fixed above:

Connectivity over Lake Lugano

All together:

Weights

# observe different location (folder with `nb` extension)
# this is dataset that 'forces' connectivity in certain areas
tg3o_nb <- 
  st_read("data/BfS/KM04-00-c-suis-2022-q/01_INST/Vegetationsfläche_vf/K4_polg20220101_vf_nb/K4polg20220101vf_ch2007Poly.shp",
          as_tibble = TRUE) %>% 
  st_transform(2056) %>% 
  rename(GMDNR = id,
         GMDNAME = name)

Reading layer `K4polg20220101vf_ch2007Poly' from data source 
  `C:\projects\ISPM_excess-mortality-voting\data\BfS\KM04-00-c-suis-2022-q\01_INST\Vegetationsfläche_vf\K4_polg20220101_vf_nb\K4polg20220101vf_ch2007Poly.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 2145 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2485203 ymin: 1075283 xmax: 2833175 ymax: 1295864
Projected CRS: CH1903+ / LV95

tg3o_wm_q <- poly2nb(tg3o_nb, queen = TRUE, row.names = tg3o$GMDNR)

tg3o_wm_q

Neighbour list object:
Number of regions: 2145 
Number of nonzero links: 11296 
Percentage nonzero weights: 0.2455105 
Average number of links: 5.2662 

summary(tg3o_wm_q)

Neighbour list object:
Number of regions: 2145 
Number of nonzero links: 11296 
Percentage nonzero weights: 0.2455105 
Average number of links: 5.2662 
Link number distribution:

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  18  19 
 48 121 234 387 452 389 233 149  73  27  13   5   3   5   1   3   1   1 
48 least connected regions:
286 298 406 412 416 590 591 598 622 631 710 757 759 820 896 911 992 1002 1115 1117 1118 1123 1138 1144 1154 1168 1171 1174 1181 1407 1499 1502 1513 1557 1580 1585 1586 1591 1903 1909 1920 1923 1950 1954 1964 2018 2070 2125 with 1 link
1 most connected region:
1702 with 19 links

write_rds(tg3o_wm_q, "data/nb/tg3o_wm_q.Rds")
nb2INLA("data/nb/tg3o_wm_q.adj", tg3o_wm_q) 

All together

Original connectivity, notice disconnected ‘islands’ of polygons:

Forcing extra connections

Clearer when only the differences are shown

List of ‘extra’, ‘forced’ connections:

• Lötschberg Tunnel connecting Kandersteg and Ferden
  
• Lötschberg Base Tunnel connecting Frutigen and Raron
  
• road connection between Ried-Brieg & Simplon

References

• Bache S, Wickham H (2022). magrittr: A Forward-Pipe Operator for R. Rpackage version 2.0.3, https://CRAN.R-project.org/package=magrittr.
• Bivand R, Nowosad J, Lovelace R (2022). spData: Datasets for SpatialAnalysis. R package version 2.2.0,https://CRAN.R-project.org/package=spData.
• Bivand R, Wong DWS (2018). “Comparing implementations of global andlocal indicators of spatial association.” TEST, 27(3), 716-748.doi:10.1007/s11749-018-0599-xhttps://doi.org/10.1007/s11749-018-0599-x.Bivand R (2022). “R Packages for Analyzing Spatial Data: A ComparativeCase Study with Areal Data.” Geographical Analysis, 54(3), 488-518.doi:10.1111/gean.12319 https://doi.org/10.1111/gean.12319.Bivand RS, Pebesma E, Gomez-Rubio V (2013). Applied spatial dataanalysis with R, Second edition. Springer, NY.https://asdar-book.org/.
• Firke S (2021). janitor: Simple Tools for Examining and Cleaning DirtyData. R package version 2.1.0,https://CRAN.R-project.org/package=janitor.
• Grolemund G, Wickham H (2011). “Dates and Times Made Easy withlubridate.” Journal of Statistical Software, 40(3), 1-25.https://www.jstatsoft.org/v40/i03/.
• Henry L, Wickham H (2022). purrr: Functional Programming Tools. Rpackage version 0.3.5, https://CRAN.R-project.org/package=purrr.
• Müller K, Wickham H (2022). tibble: Simple Data Frames. R packageversion 3.1.8, https://CRAN.R-project.org/package=tibble.
• Pebesma E (2018). “Simple Features for R: Standardized Support forSpatial Vector Data.” The R Journal, 10(1), 439-446.doi:10.32614/RJ-2018-009 https://doi.org/10.32614/RJ-2018-009,https://doi.org/10.32614/RJ-2018-009.
• Pebesma EJ, Bivand RS (2005). “Classes and methods for spatial data inR.” R News, 5(2), 9-13. https://CRAN.R-project.org/doc/Rnews/.Bivand RS, Pebesma E, Gomez-Rubio V (2013). Applied spatial dataanalysis with R, Second edition. Springer, NY.https://asdar-book.org/.
• R Core Team (2022). R: A Language and Environment for StatisticalComputing. R Foundation for Statistical Computing, Vienna, Austria.https://www.R-project.org/.
• Rinker TW, Kurkiewicz D (2018). pacman: Package Management for R.version 0.5.0, http://github.com/trinker/pacman.
• Tennekes M (2018). “tmap: Thematic Maps in R.” Journal of StatisticalSoftware, 84(6), 1-39. doi:10.18637/jss.v084.i06https://doi.org/10.18637/jss.v084.i06.
• Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis.Springer-Verlag New York. ISBN 978-3-319-24277-4,https://ggplot2.tidyverse.org.
• Wickham H (2022). forcats: Tools for Working with CategoricalVariables (Factors). R package version 0.5.2,https://CRAN.R-project.org/package=forcats.
• Wickham H (2022). stringr: Simple, Consistent Wrappers for CommonString Operations. R package version 1.4.1,https://CRAN.R-project.org/package=stringr.
• Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E,Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, TakahashiK, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to thetidyverse.” Journal of Open Source Software, 4(43), 1686.doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
• Wickham H, Bryan J (2022). readxl: Read Excel Files. R packageversion 1.4.1, https://CRAN.R-project.org/package=readxl.
• Wickham H, François R, Henry L, Müller K (2022). dplyr: A Grammar ofData Manipulation. R package version 1.0.10,https://CRAN.R-project.org/package=dplyr.
• Wickham H, Girlich M (2022). tidyr: Tidy Messy Data. R packageversion 1.2.1, https://CRAN.R-project.org/package=tidyr.
• Wickham H, Hester J, Bryan J (2022). readr: Read Rectangular TextData. R package version 2.1.3,https://CRAN.R-project.org/package=readr.
• Wickham H, Seidel D (2022). scales: Scale Functions forVisualization. R package version 1.2.1,https://CRAN.R-project.org/package=scales.
• Xie Y, Cheng J, Tan X (2022). DT: A Wrapper of the JavaScript Library’DataTables’. R package version 0.26,https://CRAN.R-project.org/package=DT.