Maintenant nous allons aborder la réalisation de cartes plus complexes. Il s’agira de réaliser des cartes d’un indicateurs sur plusieurs années, de découper par région pour avoir une meilleur visibilité. Ainsi que de réaliser des cartes de médianes et de quantiles très utilisées en climatologie. Nous allons également mieux habiller les cartes.
Nous allons utiliser quelques librairies supplémentaires en plus de celles de la pages précédentes. Tout d’abord installer le package basemap à partir de github.
library(devtools)Le chargement a nécessité le package : usethisinstall_github('Chrisjb/basemapR')Skipping install of 'basemapR' from a github remote, the SHA1 (a36d952e) has not changed since last install.
Use `force = TRUE` to force installationPuis charger les librairies (les installer si ce n’est pas le cas)
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(sf)Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUElibrary(basemapR)
library(gridExtra)
Attachement du package : 'gridExtra'
L'objet suivant est masqué depuis 'package:dplyr':
combinelibrary(ggh4x)
Attachement du package : 'ggh4x'
L'objet suivant est masqué depuis 'package:ggplot2':
guide_axis_logtickslibrary(ggpubr)
library(grid)
library(tmap)Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')Nous reprenons le même exemple que précédemment.
don<-read_delim("data/DonneeSiclima/siclima_calcul_759_20240606.csv")Rows: 1651584 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ";"
chr (2): phase, indicator
dbl (3): cell, year, raw_value
date (2): startdate, enddate
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.grille<-st_read("data/ShapefileMaille/safran.shp")Reading layer `safran' from data source
`/home/rnleroux/Documents/Travail/Outils/OutilsR/Website/ufos/data/ShapefileMaille/safran.shp'
using driver `ESRI Shapefile'
Simple feature collection with 8602 features and 23 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -4.805748 ymin: 41.36991 xmax: 9.581899 ymax: 51.08562
Geodetic CRS: WGS84(DD)Nous allons regarder comment la médiane du nombre de jour chaud en France à évoluer au cours de deux périodes, entre 1961 et 1990 (H1) et entre 1991 et 2020 (H2). Pour cela nous manipulons le jeu de données pour créer une colonne période, et ensuite calculer la médiane et les quantiles.
Cet grille est bien un objet spatial. On peut la répresenter facilement
indicatorPeriode<-don%>%
filter(indicator=="hdaystmax")%>%
mutate(periode=case_when(year>=1961 & year<=1990~"H1",
year>=1991 & year<=2020~"H2"))%>%
filter(!is.na(periode))%>%
mutate(periode=as.factor(periode))%>%
mutate(periode=fct_relevel(periode,c("H1","H2")))%>%
group_by(periode,cell)%>%
summarize(Mediane=median(raw_value),Q05=quantile(raw_value,0.05),Q95=quantile(raw_value,0.95))%>%
ungroup()`summarise()` has grouped output by 'periode'. You can override using the
`.groups` argument.indicatorPeriode# A tibble: 17,204 × 5
periode cell Mediane Q05 Q95
<fct> <dbl> <dbl> <dbl> <dbl>
1 H1 4 0 0 0
2 H1 5 0 0 0
3 H1 13 0 0 0
4 H1 14 0 0 0
5 H1 15 0 0 0
6 H1 16 0 0 0
7 H1 17 0 0 0
8 H1 18 0 0 0
9 H1 25 0 0 0
10 H1 26 0 0 0
# ℹ 17,194 more rowsNous pouvons maintenant faire la jointure avec la grille SAFRAN pour avoir la dimension spatiale
MedianeSp<-left_join(grille,indicatorPeriode)Joining with `by = join_by(cell)`Et nous pouvons passer à la figure avec encore et toujours ggplot. Nous allons rajouter un fond de carte google earth (oui c’est mal), pour avoir quelque chose de plus jolie à l’écran. Pour définir l’emprise de l’image sattelite, nous prenons les limites de notre objet spatiale.Le niveau de zoom est laissé à l’appréciation de l’utilisateur.
BoudaryBox<-st_bbox(MedianeSp) indicatorPeriodeMap<-ggplot() +
base_map(BoudaryBox, increase_zoom = 2, basemap = 'google-satellite') +
geom_sf(data = MedianeSp, mapping=aes(fill=Mediane),color=NA) +
labs(fill="Nombre de jours\n supérieurs à 35°C")+
ggtitle("Evolution du nombre de jours chauds")+
scale_fill_distiller(palette = "OrRd",direction = 1)+
facet_grid(~periode)please cite: map data © 2020 GoogleindicatorPeriodeMap
Maintenant si on veut rajouter les quantiles,et la médiane, il faut rassembler ces trois valeurs dans une seule colonne avant de faire la carte.
MedianeSpQuant<-MedianeSp%>%
pivot_longer(cols = c("Mediane","Q05","Q95"),names_to = "stat",values_to = "values")
indicatorPeriodeMapQuantile<-ggplot() +
base_map(BoudaryBox, increase_zoom = 2, basemap = 'google-satellite') +
geom_sf(data = MedianeSpQuant, mapping=aes(fill=values),color=NA) +
labs(fill="Nombre de jours\n supérieurs à 35°C")+
ggtitle("Evolution du nombre de jours chauds")+
scale_fill_distiller(palette = "OrRd",direction = 1)+
facet_grid(stat~periode)please cite: map data © 2020 GoogleindicatorPeriodeMapQuantile
Et voilà !
Maintenant nous aimerions avoir un zoom par région pour la période H2. sf et ggplot ne permettent pas de le faire simplement, nous allons donc utiliser tmap
# indicatorPeriodeMapRegion<-ggplot() +
# # base_map(BoudaryBox, increase_zoom = 2, basemap = 'google-satellite') +
# geom_sf(data = filter(MedianeSp,periode=="H2"), mapping=aes(fill=Mediane),color=NA) +
# labs(fill="Nombre de jours\n supérieurs à 35°C")+
# ggtitle("Evolution du nombre de jours chauds")+
# scale_fill_distiller(palette = "OrRd",direction = 1)+
# facet_grid(~region,scale="free")
# indicatorPeriodeMapRegion
MedianeSpRegionH2<-MedianeSp%>%
filter(periode=="H2")
tm_shape(MedianeSp) +
tm_polygons(col ="Mediane",lwd=0,
title = 'Nombre de jours chauds',
style = 'cont',
palette="OrRd")+ # continuous variable
tm_facets('region')
On complexifie encore un peu, maintenant nous allons cartographier deux indicateurs avec des unités complètements différentes. Les jours chauds avec les sommes de précipitations. On reprend notre code du départ mais cette fois on selectionne les deux indicateurs
IndicatorNames<-c("hdaystmax","rainsum")
indicatorPeriode<-don%>%
filter(indicator%in%IndicatorNames)%>%
mutate(periode=case_when(year>=1961 & year<=1990~"H1",
year>=1991 & year<=2020~"H2"))%>%
filter(!is.na(periode))%>%
mutate(periode=as.factor(periode))%>%
mutate(periode=fct_relevel(periode,c("H1","H2")))%>%
group_by(periode,cell,indicator)%>%
summarize(Mediane=median(raw_value),Q05=quantile(raw_value,0.05),Q95=quantile(raw_value,0.95))%>%
ungroup()`summarise()` has grouped output by 'periode', 'cell'. You can override using
the `.groups` argument.MedianeSp<-left_join(grille,indicatorPeriode)Joining with `by = join_by(cell)`On arrive à la partie compliquée. Comme ces indicateurs ont des unités très différentes on ne pourra pas utiliser les facet ! Mais il existe des solutions. Il faut d’abord réaliser les cartes pour ces deux indicateurs.
JoursChaudMap<-ggplot() +
base_map(BoudaryBox, increase_zoom = 2, basemap = 'google-satellite') +
geom_sf(data = filter(MedianeSp,indicator=="hdaystmax"), mapping=aes(fill=Mediane),color=NA) +
labs(fill="Nombre de jours\n supérieurs à 35°C")+
ggtitle("Evolution du nombre de jours chauds")+
scale_fill_distiller(palette = "OrRd",direction = 1)+
facet_grid(~periode)please cite: map data © 2020 GoogleJoursChaudMap
PrecipitationMap<-ggplot() +
base_map(BoudaryBox, increase_zoom = 2, basemap = 'google-satellite') +
geom_sf(data = filter(MedianeSp,indicator=="rainsum"), mapping=aes(fill=Mediane),color=NA) +
labs(fill="Précipitation (mm)")+
ggtitle("Evolution des précipitations")+
scale_fill_distiller(palette = "Blues",direction = 1)+
facet_grid(~periode)please cite: map data © 2020 GooglePrecipitationMap
Puis on les associe grace à la fonction arrangeGrob Pour informations, il existe beaucoup de packages pour assembler les graphes, cowplot, gridExtra etc…
MapSplit<-list(JoursChaudMap,PrecipitationMap)
Map<-arrangeGrob(grobs=MapSplit, ncol=1,nrow=2,common.legend = F,clip = 'inherit',x = 0, hjust =0 )
plot(Map)