Catégories
R

Minimal

Day 16 of 30DayMapChallenge : « minimal » (previously).

Bare Réunion…

# devtools::install_github("tylermorganwall/rayshader")
library(rayshader)
library(raster)
library(tidyverse)
library(fs)

# ftp://BD_ALTI_ext:docoazeecoosh1Ai@ftp3.ign.fr/BDALTIV2_2-0_25M_ASC_RGR92UTM40S-REUN89_D974_2016-03-11.7z
# We get a bunch of ASCII files : open and merge, convert to matrix used by rayshader
elmat <- dir_ls("~/BDALTIV2_2-0_25M_ASC_RGR92UTM40S-REUN89_D974_2016-03-11/BDALTIV2/1_DONNEES_LIVRAISON_2020-06-00408/BDALTIV2_MNT_25M_ASC_RGR92UTM40S_REUN89_D974/") %>% 
  map(raster) %>% 
  reduce(merge)%>% 
  raster_to_matrix()

elmat %>%
  sphere_shade(texture = "desert") %>%
  add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_3d(elmat, zscale = 10, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))

render_snapshot("~/reunion_3d.jpeg")
Bare Réunion
Catégories
R

Food

Day 15 of 30DayMapChallenge : « food » (previously).

Results of State food controls in restaurants in La Réunion.

library(dplyr)
library(sf)
library(janitor)
library(mapsf)
library(glue)
library(lubridate)

# https://geoservices.ign.fr/adminexpress 
# COG
dep <- read_sf("~/ADE-COG_2-1_SHP_WGS84G_FRA/DEPARTEMENT.shp") %>% 
  filter(INSEE_DEP == "974")

# https://dgal.opendatasoft.com/explore/dataset/export_alimconfiance/download/?format=shp
alim <- read_sf("~/../Downloads/export_alimconfiance.shp") %>% 
  clean_names() %>% 
  filter(str_sub(code_postal, 1, 3) == "974") %>% 
  mutate(date_inspec = ymd_hms(date_inspec),
         synthese_ev = factor(synthese_ev, levels = c("A corriger de manière urgente",               
                                                      "A améliorer",
                                                      "Satisfaisant",
                                                      "Très satisfaisant")))
alim %>% 
  ggplot() +
  geom_sf(data = dep,  fill = "grey90", color = "lightblue3") +
  geom_sf(aes(color = synthese_ev)) +
  scale_color_manual(values = c(
    "A corriger de manière urgente" = "red",               
    "A améliorer" = "orange",
    "Satisfaisant" = "blue",
    "Très satisfaisant" = "green")) +
  labs(title = "Food quality control",
       subtitle = "La Réunion",
       color = "Control result",
       caption = glue("data: alim-confiance.gouv.fr
                      {min(alim$date_inspec)} - {max(alim$date_inspec)}
                  r.iresmi.net - {Sys.Date()}")) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "lightblue1"),
        plot.caption = element_text(size = 7))


ggsave("food.png", width = 20, height = 20, units = "cm", scale = 1.1)
Food controls
Catégories
R

Sugar

Day 13 of 30DayMapChallenge : « 5 minutes map » (previously).

Sugar cane fields in La Réunion

library(tidyverse)
library(sf)

# RPG région La Réunion édition 2021 
# https://wxs.ign.fr/0zf5kvnyfgyss0dk5dvvq9n7/telechargement/prepackage/RPG_REGION_PACK_DIFF_2021-01-01$RPG_2-0__SHP_RGR92UTM40S_R04_2021-01-01/file/RPG_2-0__SHP_RGR92UTM40S_R04_2021-01-01.7z
parc <- read_sf("~/PARCELLES_GRAPHIQUES.shp")

# https://geoservices.ign.fr/adminexpress 
# COG
dep <- read_sf("~/DEPARTEMENT.shp") %>% 
  filter(INSEE_DEP == "974")

parc %>% 
  ggplot() +
  geom_sf(data = dep, fill = "white", color = "blue") +
  geom_sf(aes(fill = CODE_CULTU %in% c(
              "CSA",
              "CSF",
              "CSI",
              "CSP",
              "CSR")), color = NA) +
  scale_fill_manual(values = list("TRUE" = "goldenrod",
                                  "FALSE" = "lightgrey"),
                    labels = list("TRUE" = "Cane",
                                  "FALSE" = "Other")) +
  labs(title = "Sugar cane fields",
       subtitle = "La Réunion",
       fill = "Crop") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "lightblue1"))

ggsave("reunion.png", width = 20, height = 20, units = "cm", scale = 1.1)
Cane
Catégories
R

Scales

Day 12 of 30DayMapChallenge : « scales » (previously).

French « départements » at the same size but tidily ordered by area.

library(tidyverse)
library(sf)
library(units)
library(tmap)

# data from https://geoservices.ign.fr/adminexpress
fr <- read_sf("~/adminexpress.gpkg", 
              layer = "departement")

png("scales.png", width = 20, height = 12.36, units = "cm", res = 300)
fr %>%  
  mutate(surf_km2 =  set_units(st_area(geom), km^2), 
         rank = dense_rank(surf_km2), 
         nom_dep = fct_rev(fct_reorder(as_factor(nom_dep), rank))) %>% 
  tm_shape() +
  tm_polygons(fill = "lightgrey") +
  tm_facets(by = "nom_dep", free.scales = TRUE)
dev.off()  
Départements
Catégories
R

Sentinel

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.

Infrared satellite image from Termignon
Catégories
R

Get lost

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)
Get lost
Catégories
R

Use data from Wikipedia

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)

# read all tables from the page
# keep onle those containing a country column
# clean and merge
launch <- read_html("https://en.wikipedia.org/wiki/List_of_rocket_launch_sites") %>% 
  html_table() %>% 
  keep(~ names(.x)[[1]] == "Country") %>% 
  map(~ .x %>% mutate(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=robin") +
  labs(title = "Active space launch sites",
       subtitle = year(Sys.Date()),
       caption = glue("data : https://en.wikipedia.org/wiki/List_of_rocket_launch_sites
                      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("space.png", width = 20, height = 12.36, units = "cm", scale = 1.1)
Active space launch sites
Catégories
R

The giant French Olympic-size swimming pool

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 !

The giant French Olympic-size swimming pool
Catégories
R

Use data from data.gouv.fr

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)
Réseau ferré français
Catégories
R

Use data from Openstreetmap

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)

Kyiv