ASPE - Opérations géographiques
Pascal Irz
Le 09 août 2022
Objectif
Ce tutoriel montre quelques exemples de traitements géographiques à partir de la base Aspe, afin de :
- créer des objets géographiques à partir des tables de la base
- opérer des sélections géographiques (points et stations contenus dans des polygones)
- reprojeter pour homogénéiser les systèmes de coordonnées qui coexistent dans la base
- produire des cartes
Le chargement du package {aspe} est supposées déjà
réalisé (voir ce support),
de même que l’importation initiale des données.
Chargement des packages et des données
library(aspe)
library(tidyverse)
library(mapview)
load(file = "raw_data/tables_sauf_mei_2021_10_21_11_44_01.RData")Les systèmes de coordonnées (CRS)
Pour des traitements géographiques, il est nécessaire que tous les objets manipulés soient dans un même CRS.
Quelle est la fréquence de chaque CRS parmi les stations et les points de prélèvement ?
crs_stations <- station %>%
rename(typ_id = sta_typ_id) %>%
left_join(y = ref_type_projection) %>%
pull(typ_libelle_sandre) %>%
as.character() %>%
table() %>%
as.data.frame() %>%
purrr::set_names(c("CRS", "Nombre de stations"))
crs_points <- point_prelevement %>%
rename(typ_id = pop_typ_id) %>%
left_join(y = ref_type_projection) %>%
pull(typ_libelle_sandre) %>%
as.character() %>%
table() %>%
as.data.frame() %>%
purrr::set_names(c("CRS", "Nombre de points"))
crs_tot <- full_join(x = crs_points,
y = crs_stations)
DT::datatable(crs_tot, width = 600, rownames = FALSE)Reprojection des stations
La table ref_type_projection fait correspondre un code
(typ_id) à différentes manières de désigner le système de
projection (le CRS). Ici nous utiliserons en particulier le libellé
(typ_libelle_sandre) et le code EPSG de chacun
des CRS.
Concernant les CRS utilisés pour les stations et les points de
prélèvement, la table ref_type_projection contient les
données suivantes :
ref_type_projection %>%
filter(typ_libelle_sandre %in% (crs_tot %>% pull(CRS))) %>%
as.data.frame() %>%
DT::datatable()Homogénéisation des systèmes de coordonnées
Comme la base comprend les outre-mers, il faut utiliser un CRS mondial, ce qui n’est pas le cas des systèmes Lambert. On choisit donc le WGS84 qui est utilisé par les GPS, OpenStreetMaps, etc.
Ajout à la table station des codes EPSG de chacune des
stations
station <- station %>%
left_join(y = ref_type_projection,
by = c("sta_typ_id" = "typ_id"))Reprojection et extraction des coordonnées des stations en WGS84
On emploie la fonction geo_convertir_coords_df() qui
donne en sortie un dataframe avec les coordonnées
homogénéisées.
coords_wgs84 <- geo_convertir_coords_df(df = station,
var_x = sta_coordonnees_x,
var_y = sta_coordonnees_y,
var_id = sta_id,
var_crs_initial = typ_code_epsg,
crs_sortie = 4326) %>%
rename(x_wgs84 = X, y_wgs84 = Y)Les premières lignes du dataframe :
coords_wgs84 %>%
head() %>%
DT::datatable(rownames = FALSE)Ajout des coordonnées WGS84 par jointure gauche
(left_join()) et suppression des colonnes qui ne serviront
plus :
station <- station %>%
left_join(y = coords_wgs84) %>%
select(-(sta_geometrie:typ_code_epsg))Vérification que les colonnes correspondant aux coordonnées en WGS84 ont bien été créées.
names(station)## [1] "sta_id" "sta_code_sandre"
## [3] "sta_libelle_sandre" "sta_enh_id"
## [5] "sta_eligibilite_ipr_iprplus" "sta_com_code_insee"
## [7] "sta_point_km_aval" "sta_localisation_precise"
## [9] "sta_code_national_masse_eau" "x_wgs84"
## [11] "y_wgs84"A ce stade, l’objet station est un simple
dataframe avec des colonnes correspondant aux coordonnées.
Pour effectuer une sélection spatiale, il faut manipuler des objets
spatiaux comme en SIG. Ces manipulations sont relativement simples grâce
au package sf.
Création d’un objet géographique pour les stations
Il s’agit de créer l’homologue d’une couche SIG des stations avec la
fonction st_as_sf() du package sf.
station_geo <- station %>%
sf::st_as_sf(coords = c("x_wgs84", "y_wgs84"),
crs = 4326)Visualisation d’un échantillon de 1000 stations.
station_geo %>%
sample_n(1000) %>%
mapview::mapview()Il y a bien quelques points aberrants mais c’est cohérant dans l’ensemble.
Sélection géographique
Principe et données
Dans cet exemple, on sélectionne les stations situées en Bretagne. Pour réaliser l’opération, on emploie le polygone de contour de la région. Les contours des départements seront également utilisés pour produire les cartes.
On utilise le package {COGiter}
qui contient le découpage administratif du pays. Ici nous utilisons les
objets regions_metro_geo et
departements_metro_geo.
library(COGiter)
ggplot(regions_metro_geo) +
geom_sf(aes(fill = REG)) +
geom_sf(data = departements_metro_geo,
alpha = 0)
sf::st_crs(regions_metro_geo)## Coordinate Reference System:
## User input: EPSG:2154
## wkt:
## PROJCRS["RGF93 / Lambert-93",
## BASEGEOGCRS["RGF93",
## DATUM["Reseau Geodesique Francais 1993",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4171]],
## CONVERSION["Lambert-93",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",46.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",3,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",44,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",700000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",6600000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["France"],
## BBOX[41.15,-9.86,51.56,10.38]],
## ID["EPSG",2154]]Sélection
On va d’abord attribuer à chaque station son département, puis
filtrer sur les 4 départements de la région. En général, la jointure
géographique requiert que les CRS des différents objets soient
identiques, ce qui n’est pas la cas ici car les objets issus de
COGiter sont en Lambert 93, cependant la fonction
geo_attribuer() du package {aspe} gère ce cas
de figure.
# les départements sélectionnés
mes_depts <- departements_metro_geo %>%
filter(DEP %in% c("22", "29", "35", "56"))
# attribution et filtre
station_bzh <- station_geo %>%
aspe::geo_attribuer(regions_metro_geo) %>%
filter(REG == "53")Visualisation
palette <- mapviewPalette("mapviewTopoColors")
mapview::mapview(mes_depts,
zcol = "DEP",
col.regions = palette(4),
alpha.regions = 0.3,
map.types = c("OpenStreetMap",
"CartoDB.Positron",
"CartoDB.DarkMatter",
"Esri.WorldImagery",
"OpenTopoMap")) +
mapview::mapview(station_bzh,
alpha = 0.9,
cex = 0.1,
col.regions = 'blue')Reprojection des points de prélèvement
On peut procéder comme pour les stations avec d’autres tables contenant des coordonnées (longitude et latitude) et une variable de CRS.
Par exemple pour la table des points de prélèvements
point_prelevement. Comme il y a beaucoup plus de points que
de stations, c’est plus long.
# ajout au df point_prelevement du code EPSG
point_prelevement <- point_prelevement %>%
left_join(y = ref_type_projection,
by = c("pop_typ_id" = "typ_id"))
# création du df des coordonnées en XGS84
coords_wgs84 <- geo_convertir_coords_df(df = point_prelevement,
var_x = pop_coordonnees_x,
var_y = pop_coordonnees_y,
var_id = pop_id,
var_crs_initial = typ_code_epsg,
crs_sortie = 4326) %>%
rename(pop_x_wgs84 = X,
pop_y_wgs84 = Y)
# ajout au df point_prelevement des coordonnées WGS84
point_prelevement <- point_prelevement %>%
left_join(y = coords_wgs84) %>%
select(-(pop_com_code_insee_wama:pop_fog_id_cerema),
-(pop_uti_id:typ_code_epsg))
# création d'un objet "sf" des points
point_prelevement_geo <- point_prelevement %>%
sf::st_as_sf(coords = c("pop_x_wgs84", "pop_y_wgs84"),
crs = 4326)
# sélection des points qui sont dans la région
point_prelevement_bzh <- point_prelevement_geo %>%
aspe::geo_attribuer(mes_depts) %>%
filter(!is.na(DEP))
# visualisation
mapview::mapview(mes_depts,
zcol = "DEP",
col.regions = palette(4),
alpha.regions = 0.3,
map.types = c("OpenStreetMap",
"CartoDB.Positron",
"CartoDB.DarkMatter",
"Esri.WorldImagery",
"OpenTopoMap")) +
mapview::mapview(point_prelevement_bzh,
alpha = 0.9,
cex = 0.1,
col.regions = 'blue')