Étiquette : spatial

Catégories
R

Where and when (not) to eat in France ?

Results of sanitary controls in France can be found on data.gouv.fr however, only the running year is available… Thanks to @cquest@amicale.net we can access the archives since 2017.

First a global view of the dataset :

About 800 controls per week, except during the lock-down in 2020, and a slightly lower control pressure in 2021 and 2022.

Poor results (To be improved or To be corrected urgently) are stable at around 6 %, except a recent spike in 2023?

It seems that poor results increase during the year, from June to November.

This surprising periodic phenomenon is also visible by day :

So for the « when », it is : not in summer or autumn.

What about the « where »? It seems you also could be careful in some départements…

Not good in Guadeloupe, Guyane, Réunion and the southern lower Seine valley, west of Paris.

Can we see more in details ? Using a 30 km kernel smoothing :

It confirms some hot-spots west of Paris, in Alsace, in Indre, Cher, Alpes-Maritimes and between Gironde and Landes. You are safer in Paris and Bretagne…

Has it changed ?

No real trend…

Download the supporting data and R scripts :

alimconfianceTélécharger
Catégories
R

Open and merge multiple shapefiles, updated

https://flic.kr/p/yMkw9 CC by INTVGene

This old post sees a little traffic from search engines but is a mess after many editions due to the packages evolutions.

So, how can we (chose your term) append, merge, union or combine many shapefiles or other spatial vector data in 2023 with R, preferably using tidyverse functions ?

For good measure, we want to add the source file as an attribute.

First, make some data using the geopackage available here. We generate several shapefiles (one per région) in a temporary directory :

dep <- sf::read_sf("~/data/adminexpress/adminexpress_cog_simpl_000_2022.gpkg", 
                   layer = "departement")

dep |> 
  dplyr::group_by(insee_reg) |> 
  dplyr::group_walk(\(grp_data, grp_name) sf::write_sf(grp_data, 
                                                       glue::glue("~/temp/reg_{grp_name}.shp")))

This is the way

Current and concise.

fs::dir_ls("~/temp", regexp = ".*\\.shp$") |> 
  purrr::map(\(f) sf::read_sf(f) |> 
               dplyr::mutate(source = f, .before = 1)) |> 
  dplyr::bind_rows()

Approved

Recommended, as seen in the purrr help on map_dfr (that I liked better, see below) but verbose because we have to specify that we want a sf-tibble which is lost in translation.

fs::dir_ls("~/temp", regexp = ".*\\.shp$") |> 
  purrr::map(\(f) sf::read_sf(f) |> 
               dplyr::mutate(source = f, .before = 1)) |> 
  purrr::list_rbind() |>
  dplyr::as_tibble() |> 
  sf::st_sf()

Superseded

… Sadly. That’s short and understandable.

fs::dir_ls("~/temp", regexp = ".*\\.shp$") |> 
  purrr::map_dfr(\(f) sf::read_sf(f) |> 
                   dplyr::mutate(source = f, .before = 1))

Older

I liked it, too.

fs::dir_ls("~/temp", regexp = ".*\\.shp$") |> 
  dplyr::tibble(source = _) |>
  dplyr::mutate(shp = purrr::map(source, sf::read_sf)) |>
  tidyr::unnest(shp) |>
  sf::st_sf()

Oldest

If all the files share the same attributes structure.

fs::dir_ls("~/temp", regexp = ".*\\.shp$") |> 
  purrr::map(\(f) sf::read_sf(f) |> 
               dplyr::mutate(source = f, .before = 1)) |> 
  do.call(what = rbind)

or

fs::dir_ls("~/temp", regexp = ".*\\.shp$") |> 
  purrr::map(\(f) sf::read_sf(f) |> 
               dplyr::mutate(source = f, .before = 1)) |> 
  purrr::reduce(rbind)
Catégories
R

One neighbour

Today I saw this :

So, Portugal has only one terrestrial neighbour. Opencage lists the other countries with only one neighbour.

Can we get this by ourselves ?

library(tidyverse)
library(rnaturalearth)
library(spdep)

# get data from Natural Earth
countries <- ne_countries(scale = 10, returnclass = "sf") |> 
  st_make_valid() |> 
  group_by(sov_a3, sovereignt) |> 
  summarise(.groups = "drop")

# for joining codes and names later
lookup <- countries |> 
  st_drop_geometry()

# get neighbours with {spdep}
neighbours <- countries |> 
  poly2nb()

# we only keep countries with 1 neighbour
# and add names
neighbours |> 
  set_names(countries$sov_a3) |> 
  keep(\(x) length(x) == 1) |> 
  enframe("sov_a3", "nb") |> 
  unnest(nb) |> 
  filter(nb != 0) |> 
  mutate(nb = countries$sov_a3[nb]) |> 
  left_join(lookup, by = "sov_a3") |> 
  left_join(lookup, by = c("nb" = "sov_a3"), suffix = c("", "_nb")) |> 
  relocate(sovereignt, .after = 1)
# A tibble: 16 × 4
   sov_a3 sovereignt         nb    sovereignt_nb           
   <chr>  <chr>              <chr> <chr>                   
 1 BRN    Brunei             MYS   Malaysia                
 2 CAN    Canada             US1   United States of America
 3 DN1    Denmark            DEU   Germany                 
 4 DOM    Dominican Republic HTI   Haiti                   
 5 GMB    Gambia             SEN   Senegal                 
 6 HTI    Haiti              DOM   Dominican Republic      
 7 IRL    Ireland            GB1   United Kingdom          
 8 KOR    South Korea        PRK   North Korea             
 9 LSO    Lesotho            ZAF   South Africa            
10 MCO    Monaco             FR1   France                  
11 PNG    Papua New Guinea   IDN   Indonesia               
12 PRT    Portugal           ESP   Spain                   
13 QAT    Qatar              SAU   Saudi Arabia            
14 SMR    San Marino         ITA   Italy                   
15 TLS    East Timor         IDN   Indonesia               
16 VAT    Vatican            ITA   Italy

I miss Bahrain / Saudi Arabia because the Natural Earth dataset is not detailed enough…

Border on a small island

And the « new » Canada / Denmark border is also not taken into account…

While we are at it, which country has more neighbors ?

neighbours |> 
  set_names(countries$sov_a3) |> 
  unclass() |> 
  enframe("sov_a3", "nb") |> 
  unnest(nb) |> 
  count(sov_a3, sort = TRUE) |> 
  left_join(lookup, by = "sov_a3")

So, China… 

# A tibble: 209 × 3
   sov_a3     n sovereignt                      
   <chr>  <int> <chr>                           
 1 CH1       15 China                           
 2 RUS       14 Russia                          
 3 BRA       11 Brazil                          
 4 FR1       11 France                          
 5 COD        9 Democratic Republic of the Congo
 6 DEU        9 Germany                         
 7 AUT        8 Austria                         
 8 SDN        8 Sudan                           
 9 SRB        8 Republic of Serbia              
10 TUR        8 Turkey                          
# … with 199 more rows

However it includes the disputed part of « Kashmir », claimed by India and China, so it’s 14 neighbours, like Russia.

Catégories
R

Map your Strava activities

Cilaos - CC BY-NC-SA r.iresmi.net

Use the Strava API to map your outings with R using the {rStrava} package.

library(tidyverse)
library(leaflet)
library(rStrava)
library(sf)

# Get your credentials from https://www.strava.com/settings/api
app_name <- "*******"
client_id <- "*******"
client_secret <- "*******"

#' Convert Google polylines from Strava activities to {sf} polylines
#'
#' @param gp string : encoded polyline
#'
#' @return {sf} polyline
gp2sf <- function(gp) {
  gp |> 
    googlePolylines::decode() |> 
    map_dfr(
      function(df) {
        df |> 
          st_as_sf(coords = c("lon", "lat")) |> 
          st_combine() |> 
          st_cast("LINESTRING") |> 
          st_sf() 
      }) |> 
    pull(1)
}

# Get activities
activities <- httr::config(
  token = strava_oauth(
    app_name,
    client_id,
    client_secret,
    app_scope = "activity:read_all",
    cache = TRUE)) |> 
  get_activity_list() |> 
  compile_activities() 

# Map 
activities |> 
  mutate(geom = gp2sf(map.summary_polyline)) |> 
  st_sf(crs = "EPSG:4326") |> 
  leaflet() |> 
  addTiles() |> 
  addPolylines()
Catégories
R

Minimal

https://flic.kr/p/q1Ept2 - CC by Jean-Paul Tonnelier

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

https://flic.kr/p/21jc567 - CC by Frank

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

https://flic.kr/p/9ZjQMj - CC by Marcel Holyoak

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