What if all private swimming pools could be merged into one 25 m width pool? OSM is not just a map, it’s a database, so ask OSM… I know that not all swimming pools are present in OSM, but it’s just an exercise 🙂 and it can give us an order of magnitude or at least a minimum.
We will use a simplified map of France in the background and ask to the Overpass API (with {osmdata}) all « leisure=swimming_pool » and « access=private ». It takes 6 hours and 15 Go of RAM…
The CORINE Landcover dataset is distributed as a geopackage weighting more than 8 Go. To limit the memory used when we only work on a subset, we can clip it at opening time. Here we will map the Cyprus Island :
https://www.flickr.com/photos/54544400@N00 - CC by Mikael Leppä
The Suunto watches (Spartan, Suunto 9,…) can record waypoints (or POIs) but although they can be visualized in the Suunto app (or on the watch), they cannot be exported to be used with other tools. It used to be possible to access them from the Movescount website but it was discontinued a few months ago.
Now it’s just a matter of copying the file (via USB, bluetooth, etc.) to a PC, connecting to the database and finding where the POIs are stored (in a table called… pois) ; we can use QGIS or R for that…
With R, we can select this year POIs (dates are stored as UNIX timestamp) and export them as GeoJSON :
The current 2021 french administrative limits database (ADMIN EXPRESS from IGN) is more detailed than the original version (from 50 MB zipped in 2017 to 500 MB zipped now), thanks to a more detailed geometry being currently based on the BD TOPO. However we don’t always need large scale details especially for web applications. The commune layer itself is a huge 400 MB shapefile not really usable for example in a small scale leaflet map. Even the light version (ADMIN EXPRESS COG CARTO) is 120 MB for the commune layer.
Using sf::st_simplify() in R or a similar command in QGIS on these shapefiles would create holes or overlapping polygons, shapefiles not being topologically aware. We could probably convert to lines, build topology, simplify, clean, build polygons in GRASS or ArcGis, but it’s quite a hassle…
A nice solution is using Mapshaper on mapshaper.org, or better for reproducibility using {mapshaper} in R. For such large dataset it is advised to use a node.js install instead of relying on the package’s embedded version.
in red the original, in black the simplified version with départements in bold
On Debian-like :
> sudo apt-get install nodejs npm
or on windows : install https://nodejs.org/. If needed add C:\Users\xxxxxxxx\AppData\Roaming\npm to your $PATH.
> npm config set proxy "http://login:password@proxy:8080" # if necessary
> npm install -g mapshaper
For ms_simplify() we will set sys = TRUE to take advantage of the node.js executable. Experiment with the other parameters to get a resolution that suits you. Here we use Visvalingam at 3%, squeezing the commune layer from 400 MB to 30 MB. From here we rebuild departement, region and epci with ms_dissolve() commands. Then we join back with original attributes and export in a geopackage with some metadata.
library(tidyverse)
library(sf)
library(rmapshaper)
library(geojsonio)
library(janitor)
library(fs)
# ADMIN EXPRESS COG France entière édition 2021 (in WGS84)
# ftp://Admin_Express_ext:Dahnoh0eigheeFok@ftp3.ign.fr/ADMIN-EXPRESS-COG_3-0__SHP__FRA_WM_2021-05-19.7z
# also available on :
# http://files.opendatarchives.fr/professionnels.ign.fr/adminexpress/ADMIN-EXPRESS-COG_3-0__SHP__FRA_WM_2021-05-19.7z
# originals ---------------------------------------------------------------
source_ign <- "~/sig/ADMINEXPRESS/ADMIN-EXPRESS-COG_3-0__SHP__FRA_2021-05-19/ADMIN-EXPRESS-COG/1_DONNEES_LIVRAISON_2021-05-19/ADECOG_3-0_SHP_WGS84G_FRA"
com <- source_ign %>%
path("COMMUNE.shp") %>%
read_sf() %>%
clean_names()
dep <- source_ign %>%
path("DEPARTEMENT.shp") %>%
read_sf() %>%
clean_names()
reg <- source_ign %>%
path("REGION.SHP") %>%
read_sf() %>%
clean_names()
epci <- source_ign %>%
path("EPCI.shp") %>%
read_sf() %>%
clean_names()
# simplify ---------------------------------------------------------------
check_sys_mapshaper()
# 6 min
# using a conversion to geojson_json to avoid encoding problems
com_simpl <- com %>%
geojson_json(lat = "lat", lon = "long", group = "INSEE_COM", geometry = "polygon", precision = 6) %>%
ms_simplify(keep = 0.03, method = "vis", keep_shapes = TRUE, sys = TRUE)
dep_simpl <- com_simpl %>%
ms_dissolve(field = "insee_dep", sys = TRUE)
reg_simpl <- com_simpl %>%
ms_dissolve(field = "insee_reg", sys = TRUE)
epci_simpl <- com_simpl %>%
ms_dissolve(field = "siren_epci", sys = TRUE)
# add attributes and export ----------------------------------------------
destination <- "~/donnees/ign/adminexpress_simpl.gpkg"
com_simpl %>%
geojson_sf() %>%
st_write(destination, layer = "commune",
layer_options = c("IDENTIFIER=Communes Adminexpress 2021 simplifiées",
"DESCRIPTION=France WGS84 version COG (2021-05). Simplification mapshaper."))
dep_simpl %>%
geojson_sf() %>%
left_join(st_drop_geometry(dep), by = "insee_dep") %>%
st_write(destination, layer = "departement",
layer_options = c("IDENTIFIER=Départements Adminexpress 2021 simplifiés",
"DESCRIPTION=France WGS84 version COG (2021-05). Simplification mapshaper."))
reg_simpl %>%
geojson_sf() %>%
left_join(st_drop_geometry(reg), by = "insee_reg") %>%
st_write(destination, layer = "region",
layer_options = c("IDENTIFIER=Régions Adminexpress 2021 simplifiées",
"DESCRIPTION=France WGS84 version COG (2021-05). Simplification mapshaper."))
epci_simpl %>%
geojson_sf() %>%
mutate(siren_epci = str_remove(siren_epci, "200054781/")) %>% # remove Grand Paris
left_join(st_drop_geometry(epci), by = c("siren_epci" = "code_siren")) %>%
st_write(destination, layer = "epci",
layer_options = c("IDENTIFIER=EPCI Adminexpress 2021 simplifiés",
"DESCRIPTION=Établissement public de coopération intercommunale France WGS84 version COG (2021-05). Simplification mapshaper."))
Télécharger le géopackage (communes, EPCI, départements, régions simplifiés, 2021-05) en WGS84 (EPSG:4326) – 22 Mo
The {geofacet} package allows to « arrange a sequence of plots of data for different geographical entities into a grid that strives to preserve some of the original geographical orientation of the entities« .
Like the previous post, it’s interesting if you view each entity as a unit and don’t care for its real size or weight, and don’t want to spend too much time manually finding the best grid.
We will again use the sameCOVID-19 dataset. We manually add the overseas départements once we have found the right grid (by trying different seeds) and adjust Corsica position.
COVID-19 deceased in hospital, by département, for 100 000 inhab.
Hexagon tessellation using the great {geogrid} package.
The départements are the second level of administrative government in France. They neither have the same area nor the same population and this heterogeneity provides a few challenges for a fair and accurate map representation (see the post on smoothing).
However if we are just interested in the départements as units, we can use a regular grid for visualization. Since France is often called the hexagon, we could even use an hexagon tiling (a fractal map !)…
Creating the grid and conserving minimal topological relations and the general shape can be time consuming, but thanks to Geogrid it’s quite easy. The geogrid dev page provides nice examples. We will reuse our code of the COVID19 animation. The resulting GIS file is provided below.
The global shape and relations are well rendered. Deformations are quite important for the small départements around Paris, but the map remains legible.
