Day 11 of 30DayMapChallenge : « red » (previously).
Today no cartography in R but in QGIS. However the satellite image processing has been made with the {sen2R} package.
We emulate infrared photography with a combination of the 8, 4 and 3 bands from the Sentinel-2 satellites, aimed at Termignon in the Vanoise National Park.
Glaciers (and clouds) in white, grasslands in light red, forest in dark red, sparsely vegetated areas in brown reds and rocks in grey.
Day 10 of 30DayMapChallenge : « bad map » (previously).
Just resampling…
library(tidyverse)
library(sf)
library(glue)
library(ggspatial)
library(ggrepel)
pref <- read_sf("~/adminexpress/adminexpress_simpl.gpkg", layer = "commune") %>%
filter(insee_reg > "06",
str_detect(statut, "Préfecture de région|Capitale")) %>%
mutate(nom = sample(nom)) %>%
st_point_on_surface()
fr <- read_sf("~/adminexpress/adminexpress_simpl.gpkg", layer = "departement") %>%
filter(insee_reg > "06")
# map
fr %>%
ggplot() +
geom_sf(color = "grey", fill = "lightgrey") +
geom_sf(data = pref, color = "black", size = 1) +
geom_text_repel(data = pref, aes(label = nom, geometry = geom), stat = "sf_coordinates", size = 3,
bg.color = "#ffffff66",
bg.r = 0.2 ) +
coord_sf(crs = "EPSG:2154") +
labs(title = "France métropolitaine",
subtitle = "Get lost",
caption = glue("data : based on IGN Adminexpress
r.iresmi.net - {Sys.Date()}")) +
theme(plot.background = element_rect(fill = "grey10", color = NA),
panel.background = element_blank(),
panel.grid = element_line(color = "grey20"),
plot.title = element_text(color = "grey90", family = "Waree"),
plot.subtitle = element_text(color = "grey90"),
plot.caption = element_text(color = "grey90", size = 6))
ggsave("bad.png", width = 20, height = 12.36, units = "cm", scale = 1.1)
Day 9 of 30DayMapChallenge : « space » (previously).
Scrape and geolocate data from Wikipedia. We will map the active space launch sites.
library(tidyverse)
library(lubridate)
library(janitor)
library(sf)
library(glue)
library(rvest)
library(rnaturalearth)
library(showtext)
font_add_google(name = "Orbitron", family = "orbitron", regular.wt = 600)
showtext_auto()
url_data <- "https://en.wikipedia.org/wiki/List_of_rocket_launch_sites"
# read all tables from the page
# keep only those containing a country column
# clean and merge
launch <- read_html(url_data) %>%
html_table() %>%
keep(~ names(.x)[[1]] == "Country") %>%
map(~ mutate(.x, across(everything(), as.character))) %>%
reduce(bind_rows) %>%
clean_names()
# keep only active sites : roughly "YYYY-"
# and extract coordinates
active <- launch %>%
filter(str_detect(operational_date, "\\d{4}s?–(\\[.*\\])?$")) %>%
mutate(x = str_extract(coordinates, "(?<=;\\s)\\-?\\d*\\.\\d*"),
y = str_extract(coordinates, "\\-?\\d*\\.\\d*(?=;)")) %>%
select(country, location, coordinates, x, y) %>%
drop_na(x, y) %>%
st_as_sf(coords = c("x", "y"), crs = "EPSG:4326")
# map
ne_countries(scale = "small", type = "countries", returnclass = "sf") %>%
ggplot() +
geom_sf(color = "darkblue", fill = "darkblue") +
geom_sf(data = active, color = "yellow", size = 1) +
coord_sf(crs = "+proj=eqearth") +
labs(title = "Active space launch sites",
subtitle = year(Sys.Date()),
caption = glue("data : {url_data}
r.iresmi.net - {Sys.Date()}")) +
theme(text = element_text(family = "orbitron", color = "yellow", size = 20),
plot.background = element_rect(fill = "grey10", color = NA),
panel.background = element_blank(),
panel.grid = element_line(color = "grey20"),
plot.caption = element_text(size = 12, color = "grey40"))
ggsave("space.png", width = 20, height = 12.36, units = "cm", dpi = 150, scale = 1.1)
Day 8 of 30DayMapChallenge : « Openstreetmap » (previously).
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…
library(dplyr)
library(readr)
library(ggplot2)
library(sf)
library(osmdata)
library(ggspatial)
library(glue)
library(httr)
library(units)
fr_bbox <- getbb("France métropolitaine", featuretype = "area")
GET("http://r.iresmi.net/wp-content/uploads/2021/12/adminexpress_simpl.gpkg_.zip",
write_disk("~/adminexpress_simpl.gpkg_.zip"))
dir.create("~/adminexpress")
unzip("~/adminexpress_simpl.gpkg_.zip", exdir = "~/adminexpress")
fr <- read_sf("~/adminexpress/adminexpress_simpl.gpkg", layer = "region") %>%
filter(insee_reg > "06")
# 6 hours
pool <- opq(fr_bbox, timeout = 600) %>%
add_osm_feature(key = "leisure", value = "swimming_pool") %>%
add_osm_feature(key = "access", value = "private") %>%
osmdata_sf()
pool_p <- pool$osm_polygons %>%
select(osm_id)
pool_mp <- pool$osm_multipolygons %>%
select(osm_id)
rm(pool) ; gc()
# keep only pools in France (bbox was wider)
pool_fr <- pool_p %>%
bind_rows(pool_mp) %>%
st_filter(fr) %>%
mutate(area = st_area(geometry))
# what resulting dimension ?
olympic_pool_dim <- pool_fr %>%
st_drop_geometry() %>%
summarise(total_area = drop_units(sum(area, na.rm = TRUE))) %>%
mutate(width = 25,
length = total_area / width)
# rotate sf object
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
# create the pool, translate it in France and rotate
olympic_pool <- c(0, 0,
0, olympic_pool_dim$length,
olympic_pool_dim$width, olympic_pool_dim$length,
olympic_pool_dim$width, 0,
0, 0) %>%
matrix(ncol = 2, byrow = TRUE) %>%
list() %>%
st_polygon() %>%
st_sfc(crs = "EPSG:2154") *
rot(pi / 4) +
c(400000, 6300000)
# map
fr %>%
st_transform("EPSG:2154") %>%
ggplot() +
geom_sf() +
geom_sf(data = st_sf(olympic_pool, crs = "EPSG:2154"), color = "blue") +
labs(title = "The giant French Olympic-size swimming pool",
subtitle = glue("What if all private swimming pools could be merged
into one 25 m width pool?"),
caption = glue("pool data from © OpenStreetMap contributors
map IGN Adminexpress
r.iresmi.net - {Sys.Date()}")) +
annotation_scale(height = unit(0.1, "cm")) +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank())
ggsave("pool.png", width = 20, height = 12.36, units = "cm", scale = 1.1)
We get an Olympic swimming pool (25 m width) measuring at least 755 km in length !
Day 6 of 30DayMapChallenge : « network » (previously).
Using GIS data directly from data.gouv.fr : railways network of France.
library(dplyr)
library(ggplot2)
library(sf)
library(glue)
library(httr)
# https://www.data.gouv.fr/fr/datasets/fichier-de-formes-des-voies-du-reseau-ferre-national/
GET("https://www.data.gouv.fr/fr/datasets/r/71e38477-1166-4074-a532-c51f3f399b09",
write_disk("~/data.zip"))
unzip("~/data.zip")
network <- read_sf("~/fichier-de-formes-des-voies-du-reseau-ferre-national.shp") %>%
st_transform("EPSG:2154")
network %>%
ggplot() +
geom_sf() +
theme_minimal() +
theme(plot.background = element_rect(fill = "white", color = NA)) +
labs(title = "Réseau ferré français",
caption = glue("données SNCF 2020 - data.gouv.fr
r.iresmi.net - {Sys.Date()}"))
ggsave("sncf.png", width = 20, height = 12.36, units = "cm", scale = 1.1)
Day 5 of 30DayMapChallenge : « Ukraine » (previously).
Using {osmdata} to extract streets :
library(dplyr)
library(ggplot2)
library(sf)
library(osmdata)
library(ggspatial)
library(glue)
ua_bbox <- getbb("kiev, ukraine", featuretype = "city")
ua <- opq(ua_bbox) %>%
add_osm_features(features = c('"highway"="primary"',
'"highway"="secondary"',
'"highway"="tertiary"',
'"highway"="unclassified"',
'"highway"="residential"')) %>%
osmdata_sf()
ua$osm_lines %>%
bind_rows(ua$osm_multilines) %>%
st_set_crs("EPSG:4326") %>%
ggplot() +
geom_sf(color = "#FFD700", size = 0.2) +
annotation_scale(bar_cols = c("#FFD700", "#0057B7"),
line_col = "#FFD700",
text_col = "#FFD700",
height = unit(0.1, "cm")) +
coord_sf(xlim = ua_bbox[c(1, 3)],
ylim = ua_bbox[c(2, 4)]) +
labs(title = "Kyiv",
caption = glue("© OpenStreetMap contributors
r.iresmi.net - {Sys.Date()}")) +
theme_void() +
theme(plot.background = element_rect(color = NA,
fill = "#0057B7"),
plot.title = element_text(color = "#FFD700"),
plot.caption = element_text(size = 5,
color = "#FFD700"))
ggsave("kyiv.png", width = 20, height = 12.36, units = "cm", scale = 1.1)
Intersecting an area of interest with a layer at opening time
Days 3 and 4 of 30DayMapChallenge : « polygons » and « green » (previously).
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 :
library(dplyr)
library(ggplot2)
library(stringr)
library(sf)
library(rnaturalearth)
library(glue)
# Using the contour of Cyprus (enlarged) from naturalearth to clip
bb <- ne_countries(scale = "medium", country = "cyprus", returnclass = "sf") %>%
st_transform("EPSG:3035") %>%
st_buffer(90000) %>%
pull(geometry) %>%
st_as_text()
# Corine Land Cover
# download from https://land.copernicus.eu/pan-european/corine-land-cover/clc2018
# (registration required)
# passing the bounding area
cyprus_clc <- read_sf("data/U2018_CLC2018_V2020_20u1.gpkg", query = glue("
SELECT *
FROM U2018_CLC2018_V2020_20u1
WHERE st_intersects(Shape, st_polygonfromtext('{bb}'))"))
legend_colors <- list("TRUE" = "mediumaquamarine",
"FALSE" = "grey90")
legend_labs <- list("TRUE" = "forest",
"FALSE" = "other")
# exclude sea (code 523)
# and classify on forest (codes 3xx)
cyprus_clc %>%
filter(Code_18 != "523") %>%
ggplot() +
geom_sf(aes(fill = str_detect(Code_18, "^3"),
color = str_detect(Code_18, "^3"))) +
scale_fill_manual(name= "type",
values = legend_colors,
labels = legend_labs) +
scale_color_manual(name = "type",
values = legend_colors,
labels = legend_labs) +
labs(title = "Cyprus",
subtitle = "Landcover",
caption = glue("data: Copernicus CLC 2018
projection LAEA
r.iresmi.net {Sys.Date()}")) +
theme_minimal() +
theme(legend.position = "bottom",
plot.caption = element_text(size = 7))
ggsave("cyprus.png", width = 20, height = 12.36, units = "cm", scale = 1.1, type = "cairo")
library(tidyverse)
library(sf)
library(glue)
library(rnaturalearth)
library(units)
# grams of carbon dioxide-equivalents per passenger kilometer
# https://en.wikipedia.org/wiki/Fuel_economy_in_aircraft
co2_eq <- set_units(88, g/km)
# countries map from Naturalearth
countries <- ne_countries(scale = "small", returnclass = "sf")
# airport code and coordinates to geolocate itineraries
airport <- read_csv("https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat",
col_names = c("airport",
"name",
"city",
"country",
"iata",
"icao",
"latitude",
"longitude",
"altitude",
"timezone",
"dst",
"tz",
"type",
"source")) %>%
# Add Kai Tak, missig from the airport data
add_row(iata = "HKGX",
name = "Kai Tak",
city = "Hong Kong",
latitude = 22.328611,
longitude = 114.194167)
# itineraries
flight <- read_delim("from-to
LYS-LHR
LHR-LYS
LYS-BOD
LYS-BOD
LYS-BOD
LYS-BOD
BOD-LYS
BOD-LYS
BOD-LYS
LYS-BOD
BOD-LYS
BOD-LGW
LHR-JNB
CPT-JNB
JNB-LHR
LHR-ORY
BOD-ORY
CDG-HKGX
HKGX-PER
SYD-HKGX
HKGX-CDG
ORY-CAY
CAY-BEL
BEL-BSB
BSB-MAO
MAO-VVI
VVI-LPB
LPB-MAO
MAO-BEL
BEL-CAY
CAY-XAU
XAU-CAY
CAY-XAU
XAU-CAY
CAY-XAU
XAU-CAY
CAY-ORY
NCE-MXP
MXP-NCE
CDG-CAY
CAY-MPY
MPY-CAY
CAY-CDG
CDG-HKG
HKG-SYD
SYD-HKG
HKG-SYD
TLN-ORY
CDG-CPH
CPH-ORY
ORY-TLN
CDG-YYZ
YYZ-SFO
SFO-YYZ
YYZ-CDG
ORY-TLN
TLN-ORY
LYS-AMS
AMS-SHJ
SHJ-KTM
KTM-SHJ
SHJ-AMS
AMS-LYS
CDG-AUH
AUH-MCT
MCT-KTM
KTM-PKR
PKR-KTM
KTM-MCT
MCT-AUH
AUH-CDG
GVA-FCO
FCO-GVA
CDG-RUN
RUN-CDG
GVA-KEF
KEF-GVA
CDG-ARN
ARN-KRN
KRN-ARN
ARN-CDG
CDG-RUN
RUN-CDG
CDG-RUN
RUN-CDG
CDG-TLS
", delim = "-")
# geolocate
flight_geo <- flight %>%
left_join(airport, by = c("from" = "iata")) %>%
left_join(airport, by = c("to" = "iata"), suffix = c("_from", "_to"))
# create lines
flight_lines <- flight_geo %>%
mutate(line = glue("LINESTRING ({longitude_from} {latitude_from}, {longitude_to} {latitude_to})")) %>%
st_as_sf(wkt = "line", crs = "EPSG:4326")
# create great circles and compute costs
flight_geo_gc <- flight_lines %>%
st_segmentize(set_units(100, km)) %>%
mutate(distance = set_units(st_length(line), km),
co2 = set_units(distance * co2_eq, t))
# totals
total_flight <- flight_geo_gc %>%
st_drop_geometry() %>%
summarise(total_distance = sum(distance, na.rm = TRUE),
total_co2 = sum(co2, na.rm = TRUE))
# map
ggplot() +
geom_sf(data = countries, fill = "lightgrey", color = "lightgrey") +
geom_sf(data = flight_geo_gc, color = "red") +
# geom_sf(data = flight_lines, color = "blue") +
coord_sf(crs = "+proj=eqearth") +
# coord_sf(crs = "+proj=robin") +
# coord_sf(crs = "+proj=fouc") +
# coord_sf(crs = "+proj=eck1") +
# coord_sf(crs = "+proj=moll") +
# coord_sf(crs = "+proj=bonne +lat_1=10") +
# coord_sf(crs = "+proj=laea") +
labs(title = "My air travel carbon footprint 1993-2023",
subtitle = glue("{round(total_flight$total_distance, -2)} km - {round(total_flight$total_co2, 1)} teqCO₂")) +
theme_minimal()
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.
So we have to use some black magic tricks to extract them :
- Make sure Suunto app has storage permission.
- Go to settings and tap many times the version number. Logging will start.
- Go to home , pull to refresh the feed
- Wait a bit
- Go again to settings tap many times the version it will stop logging
- On your Android internal storage there will be a folder called stt
Actually, with the current app version, we get an SQLite database in :
android > data > com.stt.android.suunto > files > stt-dump > amer_app.db
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 :
library(RSQLite)
library(dplyr)
library(lubridate)
library(sf)
library(leaflet)
cnx <- dbConnect(SQLite(), "~/amer_app.db")
poi <- dbGetQuery(cnx, "SELECT * FROM pois") %>%
filter(creation >= format(ymd("2022-01-01"), "%s")) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = "EPSG:4326")
poi %>%
st_write("~/pois.geojson")
Or we can easily map them :
poi %>% leaflet() %>% addCircleMarkers() %>% addTiles()
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.
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
Script et données pour 2022 :