Dans cette application, nous souhaitons effectuer des premières requêtes géographiques sur une base PostgreSQL contenant les données du RPG. L’objectif est d’examiner les cultures à proximité d’un point géographique donné, puis de comparer la composition observée avec les compositions départementale et nationale. Dans un second temps, on verra comment développer un tableau de bord pour obtenir une visualisation interactive des résultats de ces requêtes.
Comme expliqué dans l’application précédente, la base PostgreSQL est accessible depuis tous les services interactifs du cluster. Des indications sur la manière de s’y connecter figurent sur la page. Pour rappel, le schéma public est ouvert à tous en écriture pour créer des tables temporaires. Pour ces tables, essayez d’utiliser des noms uniques pour éviter d’entrer en conflit avec d’autres équipes.
2.1 Récupération des coordonnées d’un point
Commençons par récupérer les coordonnées d’un point sur Google Maps et par les stocker dans un objet spatial grâce au package sf. Choisissez un point sur la carte et copiez ses coordonnées dans le presse-papier grâce à un clic droit. Créez ensuite deux variables lat et lon contenant la latitude et la longitude du point, ainsi qu’une variable rayon contenant le rayon souhaité pour la suite de l’analyse (en mètres).
Créez ensuite un objet point contenant les informations spatiales du point et le rayon choisi, à partir d’un data.frame et en utilisant la fonction st_as_sf. On projettera les coordonnées dans le système Lambert 93 (EPSG:2154) grâce à la fonction st_transform. Notez que les coordonnées récupérées sur Google Maps sont des coordonnées GPS, un système de projection différent.
Solution
voir le code
library(tidyverse) library(aws.s3)library(sf)library(RPostgres)library(janitor)library(kableExtra)library(leaflet)library(htmlwidgets)library(RColorBrewer)coord_gmaps <-"43.447894436406216, 1.2886163291688764"lat <-as.numeric(str_split(coord_gmaps, fixed(","), simplify =TRUE)[,1])lon <-as.numeric(str_split(coord_gmaps, fixed(","), simplify =TRUE)[,2])rayon <-10000# Création d'une table sf «point» avec les coordonnées saisies# Transformation des coordonnées en système de proj 2154 (Lambert II) point <-data.frame(lon, lat, rayon) %>%st_as_sf(coords =c("lon","lat"), crs ="EPSG:4326") %>%mutate(coord_pt_gps =st_as_text(geometry)) %>%st_transform("EPSG:2154") %>%st_sf() %>%clean_names() %>%rename(geom = geometry)
2.1.1 Connexion au serveur PostgreSQL
La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL muni de l’extension PostGIS. Connectez vous à ce serveur PostgreSQL en utilisant des variables d’environnement USER_POSTGRESQL, PASS_POSTGRESQL et HOST_POSTGRESQL bien configurées.
2.1.2 Sélection des parcelles situées autour d’un point
On souhaite maintenant requêter la base de données PostgreSQL pour récupérer les parcelles se situant dans un cercle du rayon choisi autour du point choisi. Pour cela, une possibilité est de passer par la création d’une table intermédiaire à partir de l’objet point défini précédemment (la fonction write_sf permet d’écrire des données spatiales dans une table PostGIS), puis d’utiliser la fonction ST_DWithin dans une requête. Stockez les résultats de la requête dans une variable parc_prox grâce à la fonction st_read.
Solution
voir le code
# Optionnel, suppression de la table `point` si elle existeres <-dbSendQuery(cnx, "DROP TABLE IF EXISTS public.point CASCADE;")# Ecriture de la table point dans une table PostGISwrite_sf(point, cnx, Id(schema ="public", table ="point"))# Envoi de la requête de découpage du RPG autour du point sur PostGISquery <-"SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.* FROM public.point p, rpg.parcelles r WHERE ST_DWithin(p.geom, r.geom, p.rayon);"parc_prox <-st_read(cnx, query = query)
2.1.3 Visualisation avec une carte interactive leaflet
On souhaite maintenant utiliser la librairie leaflet pour créer une visualisation interactive des données. On va vouloir afficher sur la carte interactive les libellés des cultures. Pour ceci, on récupère les groupes de cultures agrégés sur l’espace de stockage du SSP Cloud, avec un léger prétraitement, comme indiqué ici :
On va créer le widget grâce à la fonction leaflet, qui prend en argument une table sf. On peut ensuitea utiliser les fonctions addTiles (ajout d’un fond de carte) et addPolygons qui permet de personnaliser l’affichage des parcelles, notamment grâce à l’argument fillColor. Nous proposons une solution juste en dessous, mais n’hésitez pas à expérimenter !
Note
Pour utiliser leaflet, il faut que les données spatiales soient en coordonnées GPS.
Solution
voir le code
# Création d'une palette de couleurs associée au groupe de culturefactpal <-colorFactor("Paired", parc_prox$code_group)# Transformation de la projection car leaflet ne connait que le WGS 84parc_prox <- parc_prox %>%st_transform(4326)# Pour ajouter un marqueur du pointpt_mark <- point %>%st_transform(4326)# Ajout du libellé des culturesparc_prox_lib <- parc_prox %>%left_join(lib_cult %>%select(-code_groupe_culture), by =c("code_cultu"="code_culture")) # Création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la cartelabels <-sprintf("<strong>id_parcel : </strong>%s<br/> <strong>Groupe culture : </strong>%s<br/> <strong>Culture : </strong>%s<br/> <strong>Surface (ha) : </strong>%s<br/> <strong>Département : </strong>%s<br/> <strong>Commune : </strong>%s<br/>", parc_prox$id_parcel, parc_prox_lib$libelle_groupe_culture, parc_prox_lib$libelle_culture, parc_prox$surf_parc, parc_prox$insee_dep, parc_prox$nom_com) %>%lapply(htmltools::HTML)# Création de la cartecarte_parc_prox_html <-leaflet(parc_prox_lib) %>%addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%addPolygons(fillColor =~factpal(code_group),weight =2,opacity =1,color ="#ffd966",dashArray ="3",fillOpacity =0.5,highlight =highlightOptions(weight =5,color ="#A40000",dashArray ="",fillOpacity =0.0,bringToFront =TRUE),label = labels,labelOptions =labelOptions(style =list("font-weight"="normal", padding ="3px 8px"),textsize ="15px",direction ="auto",encoding="UTF-8")) %>%addMarkers(data = pt_mark, ~lon, ~lat, popup =~coord_pt_gps, label =~coord_pt_gps)# Pour sauvegarder la carte# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox.html")carte_parc_prox_html
2.1.4 Composition des parcelles agricoles récupérées
On souhaite calculer des statistiques sur les parcelles récupérées. Dans une table t1, inclure le nombre de parcelles par groupe de culture, ainsi que le nombre total de parcelles parmi l’échantillon issu de la requête. Faites la même chose sur la surface des parcelles, puis calculez la surface moyenne par parcelle et groupe de culture.
2.1.5 Comparaison avec la répartition des cultures au niveau départemental et national
On va vouloir comparer la composition des parcelles à proximité du point choisi avec la composition au niveaux départemental et national. Pour ce faire, commencer par faire une jointure spatiale sur PostGIS pour récupérer le département du point. Les géométries des départements se récupèrent avec la commande suivante :
# Couche département pour récupérer le département du pointdep <-s3read_using(FUN = sf::read_sf,layer ="departement",object ="2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",bucket ="projet-funathon",opts =list("region"="")) %>%st_transform(2154)
Pour faire la jointure spatiale, on pourra utiliser la fonction st_join.
Solution
voir le code
# Jointure spatialedf <- point %>%st_join(dep) %>%st_drop_geometry() %>%select(insee_dep)dep_pt <- df[1,1]
Pour ne pas avoir à refaire de gros calculs sur la table RPG, les statistiques départementales et nationales sont disponibles directement grâce aux commandes suivantes :
Sélectionnez les statistiques du département concerné et appariez statistiques de surfaces locales, départementales et nationales par groupe de culture dans un même data.frame à afficher, comme ci-dessous.
Solution
voir le code
# Calcul des % surfaces autour du pointstat_pt <- parc_prox %>%st_drop_geometry() %>%count(code_group, wt = surf_parc) %>%add_tally(n) %>%mutate(pct_surf_local =round(100* n / nn, 1)) %>%select(code_group, pct_surf_local) # Récupération des statistiques du département concernéstat_dep_pt <- stat_dep_pt %>%filter(insee_dep %in% dep_pt) %>%select(insee_dep, code_group, libelle_groupe_culture, pct_surf) %>%rename(pct_surf_dep = pct_surf)# Appariement des statistiques locale, départementale et nationalestat_compar <- stat_fm %>%left_join(stat_dep_pt %>%select(code_group, pct_surf_dep), by ="code_group") %>%left_join(stat_pt , by ="code_group") %>%select(libelle_groupe_culture, pct_surf_local, pct_surf_dep, pct_surf_fm) %>%arrange(desc(pct_surf_local)) %>%adorn_totals() stat_compar %>%setNames(c("Groupe de cultures","Surf. locales (%)", "Surf. départ. (%)","Surf. France m. (%)")) %>%kable(format="html",caption="<span style='font-size:medium'>Comparaison des surfaces locales, départementales et nationales</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(stat_compar), bold = T, color ="white", background ="grey")
Structure des cultures
Groupe de cultures
Surf. locales (%)
Surf. départ. (%)
Surf. France m. (%)
Blé tendre
17,6
14,3
17,7
Maïs grain et ensilage
16,9
8,6
10,0
Autres céréales
10,6
12,4
3,9
Gel (surfaces gelées sans production)
8,7
4,4
1,6
Tournesol
8,3
12,6
2,5
Prairies permanentes
7,3
17,9
27,7
Autres oléagineux
7,2
3,1
0,7
Orge
5,5
3,5
6,2
Fourrage
4,7
5,9
3,7
Prairies temporaires
4,6
4,2
5,3
Protéagineux
3,1
1,8
1,2
Colza
2,7
1,8
3,5
Divers
1,1
1,3
1,3
Légumineuses à grains
0,9
0,3
0,2
Légumes ou fleurs
0,3
0,3
1,5
Estives et landes
0,2
6,9
8,1
Vignes
0,2
0,4
2,2
Autres cultures industrielles
0,2
0,1
1,8
Vergers
0,0
0,1
0,4
Plantes à fibres
NA
0,0
0,4
Riz
NA
NA
0,0
Fruits à coque
NA
0,0
0,2
Oliviers
NA
0,0
0,0
Total
100,1
99,9
100,1
2.1.6 Graphique de comparaison des cultures au niveau local, départemental et national
On souhaite faire un graphique avec ggplot2 pour visualiser la comparaison établie ci-dessus. Utilisez la fonction geom_col pour obtenir le graphique affiché ci-dessous.
Solution
voir le code
# Sélection des 10 groupes de cultures les plus répandus au niveau local tab <- stat_compar %>%filter(libelle_groupe_culture !="Total") %>%slice_head(n=10) %>%rename(local = pct_surf_local, departement = pct_surf_dep, france = pct_surf_fm)# Transposition de la table pour rassembler toutes les valeurs dans une seule variable valuetab_piv <- tab %>%pivot_longer(!libelle_groupe_culture) %>%rename(secteur = name) # Valeurs manquantes àtab_piv[is.na(tab_piv)] <-0# On réordonne les secteurs dans le "bon" ordre, avec factortab_piv$secteur <-factor( tab_piv$secteur,levels =c("france", "departement", "local"))tab_piv <- tab_piv %>%arrange(desc(secteur), desc(value))# On réordonne les cultures par surface décroissante au niveau local, avec factorx <- tab_piv %>%filter(secteur =="local") %>%arrange(value) %>%select(libelle_groupe_culture)y <-pull(x, libelle_groupe_culture)tab_piv$libelle_groupe_culture <-factor(tab_piv$libelle_groupe_culture, levels = y)# Visualisation avec `geom_col`p <-ggplot(tab_piv, aes(x = libelle_groupe_culture,y = value, fill =factor( secteur,levels =c("france", "departement", "local")))) +geom_col(position ="dodge") +labs(title ="Surfaces comparées des 10\nprincipales cultures locales, en %", x="Culture", y ="%", fill ="Secteur") +theme_classic()# Flip du graphique pour avoir des barres horizontales p +coord_flip()
2.1.7 Graphique par secteur
La fonction facet_wrap permet d’afficher plusieurs graphiques côte-à-côte. En utilisant cette fonction, réalisez le graphique affiché ci-dessous.
Solution
voir le code
# Visualisation avec `geom_col` et `facet_wrap` ggplot(tab_piv, aes(x = libelle_groupe_culture,y = value)) +geom_col(fill ="lightblue", colour ="black", position ="dodge") +labs(title ="Surface par culture", x="Culture", y ="%", fill ="Secteur") +geom_text(aes(label = value), hjust =-0.3, size =8/.pt, colour ="black") +theme_classic() +coord_flip() +facet_wrap(~secteur, nrow=3, scales='free')
2.2 Création d’un dashboard de visualisation
On souhaiterait intégrer les analyses faites ci-dessus à un tableau de bord qui offrirait des visualisations interactives. Une solution est d’utiliser shiny, une librairie qui permet la création de telles applications avec R notamment.
Une application shiny simple a deux composants principaux :
Un objet d’interface utilisateur (UI) qui contrôle la disposition et l’apparence du tableau de bord ;
Une fonction serveur qui contient les instructions nécessaires au fonctionnement de l’application.
Lorsque ces deux composants sont définis, l’application est lancée à l’aide d’une simple instruction runApp. Ici, on aimerait bien aller vers une première version de tableau de bord, qui afficherait au départ une carte interactive et laisserait aussi à l’utilisateur de choisir un rayon. Une fois ce rayon choisi et un point sélectionné grâce à un clic sur la carte, on souhaiterait (comme affiché ci-dessous) :
Afficher sur la carte les parcelles se situant à une distance inférieure au rayon choisi du point spécifié ;
Afficher des statistiques sur les parcelles en questions.
Ainsi notre UI sera un objet fluidPage composé de trois éléments, et dans la fonction de serveur il y aura deux évènements distincts à observer. Pour développer l’application, créez un répertoire my_app et à l’intérieur 3 fichiers, un fichier ui.R, un fichier server.R et un fichier utils.R (dans lequel se trouveront les fonctions utilitaires).
2.3 UI
Le fichier UI aura la forme suivante :
library(shiny)# Define UIui <-fluidPage( ...)
A l’intérieur de la fonction fluidPage et grâce à la documentation de Shiny, ajoutez les 3 éléments qui constitueront le tableau de bord : un output de type “carte leaflet”, un champ d’input numérique, et un output de type “table”.
Solution
voir le code
library(shiny)# Define UIui <-fluidPage(leafletOutput("map", height =800),numericInput("buffer_radius", "Rayon (en km) :", value =5),tableOutput("table"))
2.4 Serveur
Le code de la fonction serveur figure ci-dessous dans son intégralité, et l’objectif de cette partie va être d’implémenter dans utils.R les différentes fonctions appelées dans le code :
connect_to_db : une fonction qui renvoie une connexion à la base de données ;
query_db : une fonction qui interroge la base de données pour récupérer les parcelles se situant à l’intérieur d’un cercle de rayon radius d’un point défini par une latitude et une longitude ;
plot_surroundings : une fonction qui prend en entrée la sortie de la fonction query_db et une carte leaflet et qui ajoute une couche comportant les polygones des parcelles concernées ;
compute_stats : une fonction qui prend en entrée la sortie de la fonction query_db et renvoie des statistiques sur les parcelles concernées.
library(shiny)library(leaflet)# Définition du serveurserver <-function(input, output) {# Rendre la carte output$map <-renderLeaflet({leaflet() %>%addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%setView(lng =-1.4932577154775046, lat =46.46946179131805, zoom =12) })# Connexion à la base de données cnx <-connect_to_db()# Initialisation d'une "reactive value" pour le point sélectionné selectedPoint <-reactiveValues(lat =NULL, lng =NULL)# Gestion de l'évènement "clic"observeEvent(input$map_click, { clickData <- input$map_clickif (!is.null(clickData)) {# Stockage des coordonnées du point selectedPoint$lat <- clickData$lat selectedPoint$lng <- clickData$lng buffer_radius <- input$buffer_radius sf <-query_db(cnx, selectedPoint$lat, selectedPoint$lng, buffer_radius)# Mise à jour de la carteleafletProxy("map") %>%clearMarkers() %>%clearShapes() %>%addMarkers(lng = selectedPoint$lng, lat = selectedPoint$lat) %>%plot_surroundings(sf)# Calculs sur les données de parcelles df <-compute_stats(sf)# Update de la table affichée output$table <-renderTable({ df }) } })observeEvent(input$buffer_radius, {# Vérification qu'un point a été sélectionnéif (!is.null(selectedPoint$lat) &&!is.null(selectedPoint$lng)) {# Requête avec le nouveau rayon buffer_radius <- input$buffer_radius sf <-query_db(cnx, selectedPoint$lat, selectedPoint$lng, buffer_radius)# Mise à jour de la carteleafletProxy("map") %>%clearShapes() %>%plot_surroundings(sf)# Calculs sur les données de parcelles df <-compute_stats(sf)# Update de la table affichée output$table <-renderTable({ df }) } })}
Proposer une implémentation des fonctions connect_to_db, query_db et compute_stats décrites ci-dessus et utilisées dans la fonction “serveur”.
2.4.1connect_to_db
La fonction connect_to_db renvoie une connexion à la base de données PostgreSQL, comme vu précédemment. Les identifiants (user et password sont respectivement stockés dans les variables d’environnement USER_POSTGRESQL et PASS_POSTGRESQL.)
Solution
voir le code
#' Connection au serveur PostgreSQL. Le mot de passe doit être stocké dans la#' variable d'environnement PASS_POSTGRESQL.#' #' @returns Connexion au serveur.connect_to_db <-function() {# Connection à PostgreSQL cnx <-dbConnect(Postgres(),user =Sys.getenv("USER_POSTGRESQL"),password =Sys.getenv("PASS_POSTGRESQL"),host ="postgresql-100400.projet-funathon",dbname ="defaultdb",port =5432,check_interrupts =TRUE)return(cnx)}
2.4.2query_db
On souhaite implémenter la fonction suivante :
#' Requête la table `parcelles` pour récupérer les parcelles qui se situent#' dans un certain rayon autour d'un point repéré par une latitude et #' une longitude.#' #' @param cnx Connexion à PostgreSQL.#' @param lat Latitude.#' @param lng Longitude.#' @param radius Rayon.#' @returns Objet `sf` avec les données des parcelles concernées.query_db <-function(cnx, lat, lng, radius) { ...}
On pourra notamment utiliser plusieurs fonctions PostGIS :
ST_MakePoint qui permet de créer une géométrie POINT;
ST_SetSRID qui permet de définir le système de coordonnées pour une géométrie;
ST_Buffer qui calcule un POLYGON ou un MULTIPOLYGON qui représente tous les points dont la distance par rapport à une géométrie/géographie est inférieure ou égale à une distance donnée;
ST_Intersects qui compare deux géométries et renvoie true si elles ont une intersection non-nulle;
ainsi que des fonction de la librairie sf.
Solution
voir le code
#' Requête la table `parcelles` pour récupérer les parcelles qui se situent#' dans un certain rayon autour d'un point repéré par une latitude et #' une longitude.#' #' @param cnx Connexion à PostgreSQL.#' @param lat Latitude.#' @param lng Longitude.#' @param radius Rayon.#' @returns Objet `sf` avec les données des parcelles concernées.query_db <-function(cnx, lat, lng, radius) {# Les données spatiales sur PostgreSQL sont stockées en Lambert 93.# Pour faire le join on veut donc projeter les coordonnées `lat`` et `lng` postgis_crs <-"EPSG:2154" coordinates <-data.frame(lng =c(lng), lat =c(lat)) %>%st_as_sf(coords =c("lng", "lat"), remove =TRUE) %>%st_set_crs("EPSG:4326") %>%st_transform(postgis_crs)# Requête PostgreSQL query <-sprintf("SELECT * FROM rpg.parcelles WHERE ST_Intersects(geom, ST_Buffer(ST_SetSRID(ST_MakePoint(%f, %f), 2154), %.0f));",st_coordinates(coordinates)[1],st_coordinates(coordinates)[2], radius*1000)# Récupération des résultats sf <-st_read( cnx,query = query )return(sf)}
2.4.3compute_stats
On souhaite implémenter la fonction suivante :
#' Crée la table à afficher sur l'application grâce à des calculs sur les#' données requêtées depuis PostgreSQL.#' #' @param sf Données spatiales.#' @returns data.frame à afficher.compute_stats <-function(sf) { ...}
qui calcule sur une table spatiale renvoyée par la fonction query_db (ensemble des parcelles se situant autour d’un point donné) les statistiques suivantes par groupe de culture agregé :
Nombre de parcelles;
Pourcentage des parcelles parmi toutes les parcelles;
Surface des parcelles;
Surface des parcelles rapportée à la surface de toutes les parcelles;
Surface moyenne d’une parcelle.
On récupère les groupes de cultures agrégés sur l’espace de stockage du SSP Cloud et avec un léger prétraitement, comme précédemment :
#' Crée la table à afficher sur l'application grâce à des calculs sur les#' données requêtées depuis PostgreSQL.#' #' @param sf Données spatiales.#' @returns data.frame à afficher.compute_stats <-function(sf) { df <- sf %>%st_drop_geometry() %>%count(code_group, name ="parcelles_grp") %>%add_tally(parcelles_grp, name ="parcelles_tot") %>%mutate(pct_parcelles =round(100*parcelles_grp/parcelles_tot, 1)) %>%select(-parcelles_tot) %>%cbind(# Comptage des surfaces sf %>%st_drop_geometry() %>%count(code_group, wt = surf_parc, name ="surface_grp") %>%add_tally(surface_grp, name ="surface_tot") %>%mutate(surface_pct =round(100*surface_grp/surface_tot, 1)) %>%select(-surface_tot) %>%select(surface_grp, surface_pct) ) %>%left_join(lib_group_cult, by =c("code_group"="code_groupe_culture")) %>%select(code_group, libelle_groupe_culture, everything()) %>%arrange(desc(surface_grp)) %>%adorn_totals() %>%mutate(mean_surface =round(surface_grp/parcelles_grp, 1))return( df %>%select(-code_group) %>%setNames(c("Groupe de cultures","Nombre de parcelles","Pourcentage de parcelles","Surface (ha)","Surface (%)","Surface moyenne d'une parcelle (ha)")) )}
2.4.4plot_surroundings
La fonction plot_surroundings prend en entrée une carte leaflet ainsi qu’un objet sf récupéré grâce à query_db et renvoie une carte où on a rajouté une couche avec les polygones correspondant aux parcelles. Elle est donné ci-dessous. Essayez de comprendre à quoi correspondent les différents arguments de la fonction addPolygons.
# Création d'une palette de couleurs associée au groupe de culturepal <-brewer.pal(12, "Paired")pal <-colorRampPalette(pal)(24)factpal <-colorFactor(pal, lib_group_cult$code_groupe_culture)#' Rajoute les données d'un objet `sf` sous forme de polygones à une#' carte `leaflet`.#' #' @param leaflet_proxy Carte.#' @param sf Données spatiales.#' @returns Carte enrichie.plot_surroundings <-function(leaflet_proxy, sf) {# Transformation de la projection (WGS 84) sf <- sf %>%st_transform(4326)# Ajout des libellés des cultures sf <- sf %>%left_join(lib_cult %>%select(-code_groupe_culture), by =c("code_cultu"="code_culture")) # Création des labels à afficher au passage de la souris sur la carte. labels <-sprintf("<strong>Identifiant de la parcelle : </strong>%s<br/> <strong>Groupe culture : </strong>%s<br/> <strong>Culture : </strong>%s<br/> <strong>Surface (ha) : </strong>%s<br/> <strong>Département : </strong>%s<br/> <strong>Commune : </strong>%s<br/>", sf$id_parcel, sf$libelle_groupe_culture, sf$libelle_culture, sf$surf_parc, sf$insee_dep, sf$nom_com) %>%lapply(htmltools::HTML)return( leaflet_proxy %>%addPolygons(data = sf,fillColor =~factpal(code_group),weight =2,opacity =1,color ="#ffd966",dashArray ="3",fillOpacity =0.5,highlight =highlightOptions(weight =5,color ="#A40000",dashArray ="",fillOpacity =0.0,bringToFront =TRUE),label = labels,labelOptions =labelOptions(style =list("font-weight"="normal", padding ="3px 8px"),textsize ="15px",direction ="auto",encoding="UTF-8")) )}
2.5 Lancement de l’application
L’application RShiny se lance localement grâce à la commande shiny::runApp("my-app"). Vérifiez que tout fonctionne bien ! Une application fonctionnelle se trouve dans le répertoire app si besoin.
Code source
---title: "Première manipulation du RPG"pagetitle: "Première manipulation du RPG"keywords: ["RPG", "PostgreSQL", "PostGIS", "manipulation"]number-sections: trueknitr: opts_chunk: dev: "ragg_png" out.width: 100%---Dans cette application, nous souhaitons effectuer des premières requêtes géographiques sur une base PostgreSQL contenant les données du RPG. L'objectif est d'examiner les cultures à proximité d'un point géographique donné, puis de comparer la composition observée avec les compositions départementale et nationale. Dans un second temps, on verra comment développer un tableau de bord pour obtenir une visualisation interactive des résultats de ces requêtes. Comme expliqué dans l'application précédente, la base PostgreSQL est accessible depuis tous les services interactifs du cluster. Des indications sur la manière de s'y connecter figurent sur la page. Pour rappel, le schéma `public` est ouvert à tous en écriture pour créer des tables temporaires. Pour ces tables, essayez **d'utiliser des noms uniques** pour éviter d'entrer en conflit avec d'autres équipes.## Récupération des coordonnées d'un pointCommençons par récupérer les coordonnées d'un point sur Google Maps et par les stocker dans un objet spatial grâce au package [`sf`](https://r-spatial.github.io/sf/). Choisissez un point sur [la carte](https://www.google.fr/maps) et copiez ses coordonnées dans le presse-papier grâce à un clic droit. Créez ensuite deux variables `lat` et `lon` contenant la latitude et la longitude du point, ainsi qu'une variable `rayon` contenant le rayon souhaité pour la suite de l'analyse (en mètres).Créez ensuite un objet `point` contenant les informations spatiales du point et le rayon choisi, à partir d'un `data.frame` et en utilisant la fonction [`st_as_sf`](https://r-spatial.github.io/sf/reference/st_as_sf.html). On projettera les coordonnées dans le système Lambert 93 (`EPSG:2154`) grâce à la fonction `st_transform`. Notez que les coordonnées récupérées sur Google Maps sont des coordonnées GPS, un système de projection différent.::: {.callout-note icon=false}## Solution```{r}#| label: choix coordonnées + rayonlibrary(tidyverse) library(aws.s3)library(sf)library(RPostgres)library(janitor)library(kableExtra)library(leaflet)library(htmlwidgets)library(RColorBrewer)coord_gmaps <-"43.447894436406216, 1.2886163291688764"lat <-as.numeric(str_split(coord_gmaps, fixed(","), simplify =TRUE)[,1])lon <-as.numeric(str_split(coord_gmaps, fixed(","), simplify =TRUE)[,2])rayon <-10000# Création d'une table sf «point» avec les coordonnées saisies# Transformation des coordonnées en système de proj 2154 (Lambert II) point <-data.frame(lon, lat, rayon) %>%st_as_sf(coords =c("lon","lat"), crs ="EPSG:4326") %>%mutate(coord_pt_gps =st_as_text(geometry)) %>%st_transform("EPSG:2154") %>%st_sf() %>%clean_names() %>%rename(geom = geometry)```:::### Connexion au serveur PostgreSQLLa table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL muni de l'extension PostGIS. Connectez vous à ce serveur PostgreSQL en utilisant des variables d'environnement `USER_POSTGRESQL`, `PASS_POSTGRESQL` et `HOST_POSTGRESQL` bien configurées.::: {.callout-note icon=false}## Solution```{r}#| label: connection PostgreSQL# Connexion à PostgreSQLcnx <-dbConnect(Postgres(),user =Sys.getenv("USER_POSTGRESQL"),password =Sys.getenv("PASS_POSTGRESQL"),host =Sys.getenv("HOST_POSTGRESQL"),dbname ="defaultdb",port =5432,check_interrupts =TRUE)```:::### Sélection des parcelles situées autour d'un pointOn souhaite maintenant requêter la base de données PostgreSQL pour récupérer les parcelles se situant dans un cercle du rayon choisi autour du point choisi. Pour cela, une possibilité est de passer par la création d'une table intermédiaire à partir de l'objet `point` défini précédemment (la fonction `write_sf` permet d'écrire des données spatiales dans une table `PostGIS`), puis d'utiliser la fonction [`ST_DWithin`](https://postgis.net/docs/ST_DWithin.html) dans une requête. Stockez les résultats de la requête dans une variable `parc_prox` grâce à la fonction `st_read`.::: {.callout-note icon=false}## Solution```{r}#| label: requete SQL selection parcelles#| warnings : false# Optionnel, suppression de la table `point` si elle existeres <-dbSendQuery(cnx, "DROP TABLE IF EXISTS public.point CASCADE;")# Ecriture de la table point dans une table PostGISwrite_sf(point, cnx, Id(schema ="public", table ="point"))# Envoi de la requête de découpage du RPG autour du point sur PostGISquery <-"SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.* FROM public.point p, rpg.parcelles r WHERE ST_DWithin(p.geom, r.geom, p.rayon);"parc_prox <-st_read(cnx, query = query)```:::### Visualisation avec une carte interactive leafletOn souhaite maintenant utiliser la librairie `leaflet` pour créer une visualisation interactive des données. On va vouloir afficher sur la carte interactive les libellés des cultures. Pour ceci, on récupère les groupes de cultures agrégés sur l'espace de stockage du SSP Cloud, avec un léger prétraitement, comme indiqué ici :```{r}#| label: lib-cult-first#| code-fold: false# Récupération des libellés des différentes cultureslib_cult <-s3read_using(FUN = read_csv2,object ="2023/sujet2/diffusion/ign/rpg/REF_CULTURES_GROUPES_CULTURES_2020.csv",col_types =cols(.default =col_character()),bucket ="projet-funathon",opts =list("region"="")) %>%clean_names()lib_group_cult <- lib_cult %>%select(code_groupe_culture, libelle_groupe_culture) %>%distinct(code_groupe_culture, .keep_all=T)lib_group_cult %>%kable()```On va créer le widget grâce à la fonction [`leaflet`](https://rstudio.github.io/leaflet/), qui prend en argument une table `sf`. On peut ensuitea utiliser les fonctions [`addTiles`](https://rstudio.github.io/leaflet/basemaps.html) (ajout d'un fond de carte) et `addPolygons` qui permet de personnaliser l'affichage des parcelles, notamment grâce à l'argument `fillColor`. Nous proposons une solution juste en dessous, mais n'hésitez pas à expérimenter ! ::: {.callout-note icon=false}## NotePour utiliser `leaflet`, il faut que les données spatiales soient en coordonnées GPS.:::::: {.callout-note icon=false}## Solution```{r}#| label: affichage-parcelles# Création d'une palette de couleurs associée au groupe de culturefactpal <-colorFactor("Paired", parc_prox$code_group)# Transformation de la projection car leaflet ne connait que le WGS 84parc_prox <- parc_prox %>%st_transform(4326)# Pour ajouter un marqueur du pointpt_mark <- point %>%st_transform(4326)# Ajout du libellé des culturesparc_prox_lib <- parc_prox %>%left_join(lib_cult %>%select(-code_groupe_culture), by =c("code_cultu"="code_culture")) # Création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la cartelabels <-sprintf("<strong>id_parcel : </strong>%s<br/> <strong>Groupe culture : </strong>%s<br/> <strong>Culture : </strong>%s<br/> <strong>Surface (ha) : </strong>%s<br/> <strong>Département : </strong>%s<br/> <strong>Commune : </strong>%s<br/>", parc_prox$id_parcel, parc_prox_lib$libelle_groupe_culture, parc_prox_lib$libelle_culture, parc_prox$surf_parc, parc_prox$insee_dep, parc_prox$nom_com) %>%lapply(htmltools::HTML)# Création de la cartecarte_parc_prox_html <-leaflet(parc_prox_lib) %>%addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%addPolygons(fillColor =~factpal(code_group),weight =2,opacity =1,color ="#ffd966",dashArray ="3",fillOpacity =0.5,highlight =highlightOptions(weight =5,color ="#A40000",dashArray ="",fillOpacity =0.0,bringToFront =TRUE),label = labels,labelOptions =labelOptions(style =list("font-weight"="normal", padding ="3px 8px"),textsize ="15px",direction ="auto",encoding="UTF-8")) %>%addMarkers(data = pt_mark, ~lon, ~lat, popup =~coord_pt_gps, label =~coord_pt_gps)# Pour sauvegarder la carte# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox.html")carte_parc_prox_html```:::### Composition des parcelles agricoles récupéréesOn souhaite calculer des statistiques sur les parcelles récupérées. Dans une table `t1`, inclure le nombre de parcelles par groupe de culture, ainsi que le nombre total de parcelles parmi l'échantillon issu de la requête. Faites la même chose sur la surface des parcelles, puis calculez la surface moyenne par parcelle et groupe de culture.::: {.callout-note icon=false}## Solution```{r}#| label: stats sur les groupes de cultures #| tbl-cap: Structure des cultures au niveau local t1 <- parc_prox %>%st_drop_geometry() %>%count(code_group) %>%add_tally(n) %>%mutate(n_pct =round(100* n / nn, 1)) %>%select(-nn) %>%rename(n_parcelles = n) %>%cbind(# Surfaces parc_prox %>%st_drop_geometry() %>%count(code_group, wt = surf_parc) %>%add_tally(n) %>%mutate(surf_pct =round(100* n / nn, 1)) %>%select(-nn) %>%rename(surf_parc_ha = n) %>%select(surf_parc_ha, surf_pct) ) %>%left_join(lib_group_cult, by =c("code_group"="code_groupe_culture")) %>%select(code_group, libelle_groupe_culture, everything()) %>%arrange(desc(surf_parc_ha)) %>%adorn_totals() %>%mutate(taille_moy_parc =round(surf_parc_ha / n_parcelles, 1))t1 %>%setNames(c("Code", "Groupe de cultures", "Nombre de parcelles", "(%)", "Surface (ha)", "Surface (%)", "Taille moyenne (ha)")) %>%kable(format="html",caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(t1), bold = T, color ="white", background ="grey")```:::### Comparaison avec la répartition des cultures au niveau départemental et nationalOn va vouloir comparer la composition des parcelles à proximité du point choisi avec la composition au niveaux départemental et national. Pour ce faire, commencer par faire une jointure spatiale sur `PostGIS` pour récupérer le département du point. Les géométries des départements se récupèrent avec la commande suivante :```{r}#| label: dep-layer#| code-fold: false# Couche département pour récupérer le département du pointdep <-s3read_using(FUN = sf::read_sf,layer ="departement",object ="2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",bucket ="projet-funathon",opts =list("region"="")) %>%st_transform(2154)```Pour faire la jointure spatiale, on pourra utiliser la fonction `st_join`.::: {.callout-note icon=false}## Solution```{r}#| label: sf-spatial-join# Jointure spatialedf <- point %>%st_join(dep) %>%st_drop_geometry() %>%select(insee_dep)dep_pt <- df[1,1]```:::Pour ne pas avoir à refaire de gros calculs sur la table RPG, les statistiques départementales et nationales sont disponibles directement grâce aux commandes suivantes :```{r}#| label: stats-agg#| code-fold: false# Récupération des statistiques départementalesstat_dep_pt <-s3read_using(FUN = readr::read_rds,object ="2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",bucket ="projet-funathon",opts =list("region"=""))# Récupération des statistiques nationalesstat_fm <-s3read_using(FUN = readr::read_csv,object ="2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",col_types =cols(code_group =col_character()),bucket ="projet-funathon",opts =list("region"="")) %>%select(code_group, libelle_groupe_culture, pct_surf) %>%rename(pct_surf_fm = pct_surf)```Sélectionnez les statistiques du département concerné et appariez statistiques de surfaces locales, départementales et nationales par groupe de culture dans un même `data.frame` à afficher, comme ci-dessous.::: {.callout-note icon=false}## Solution```{r}#| label: compare-structures#| tbl-cap: Structure des cultures # Calcul des % surfaces autour du pointstat_pt <- parc_prox %>%st_drop_geometry() %>%count(code_group, wt = surf_parc) %>%add_tally(n) %>%mutate(pct_surf_local =round(100* n / nn, 1)) %>%select(code_group, pct_surf_local) # Récupération des statistiques du département concernéstat_dep_pt <- stat_dep_pt %>%filter(insee_dep %in% dep_pt) %>%select(insee_dep, code_group, libelle_groupe_culture, pct_surf) %>%rename(pct_surf_dep = pct_surf)# Appariement des statistiques locale, départementale et nationalestat_compar <- stat_fm %>%left_join(stat_dep_pt %>%select(code_group, pct_surf_dep), by ="code_group") %>%left_join(stat_pt , by ="code_group") %>%select(libelle_groupe_culture, pct_surf_local, pct_surf_dep, pct_surf_fm) %>%arrange(desc(pct_surf_local)) %>%adorn_totals() stat_compar %>%setNames(c("Groupe de cultures","Surf. locales (%)", "Surf. départ. (%)","Surf. France m. (%)")) %>%kable(format="html",caption="<span style='font-size:medium'>Comparaison des surfaces locales, départementales et nationales</span>",format.args =list(decimal.mark =",", big.mark =" "),booktabs =TRUE) %>%kable_styling(font_size =15) %>%gsub("font-size: initial !important;","font-size: 20pt !important;",.)%>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(nrow(stat_compar), bold = T, color ="white", background ="grey")```:::### Graphique de comparaison des cultures au niveau local, départemental et nationalOn souhaite faire un graphique avec `ggplot2` pour visualiser la comparaison établie ci-dessus. Utilisez la fonction `geom_col` pour obtenir le graphique affiché ci-dessous.::: {.callout-note icon=false}## Solution```{r}#| label: faire un graphique comparant les structures# Sélection des 10 groupes de cultures les plus répandus au niveau local tab <- stat_compar %>%filter(libelle_groupe_culture !="Total") %>%slice_head(n=10) %>%rename(local = pct_surf_local, departement = pct_surf_dep, france = pct_surf_fm)# Transposition de la table pour rassembler toutes les valeurs dans une seule variable valuetab_piv <- tab %>%pivot_longer(!libelle_groupe_culture) %>%rename(secteur = name) # Valeurs manquantes àtab_piv[is.na(tab_piv)] <-0# On réordonne les secteurs dans le "bon" ordre, avec factortab_piv$secteur <-factor( tab_piv$secteur,levels =c("france", "departement", "local"))tab_piv <- tab_piv %>%arrange(desc(secteur), desc(value))# On réordonne les cultures par surface décroissante au niveau local, avec factorx <- tab_piv %>%filter(secteur =="local") %>%arrange(value) %>%select(libelle_groupe_culture)y <-pull(x, libelle_groupe_culture)tab_piv$libelle_groupe_culture <-factor(tab_piv$libelle_groupe_culture, levels = y)# Visualisation avec `geom_col`p <-ggplot(tab_piv, aes(x = libelle_groupe_culture,y = value, fill =factor( secteur,levels =c("france", "departement", "local")))) +geom_col(position ="dodge") +labs(title ="Surfaces comparées des 10\nprincipales cultures locales, en %", x="Culture", y ="%", fill ="Secteur") +theme_classic()# Flip du graphique pour avoir des barres horizontales p +coord_flip()```:::### Graphique par secteurLa fonction [`facet_wrap`](https://ggplot2.tidyverse.org/reference/facet_wrap.html) permet d'afficher plusieurs graphiques côte-à-côte. En utilisant cette fonction, réalisez le graphique affiché ci-dessous.::: {.callout-note icon=false}## Solution```{r}#| label: facet-wrap # Visualisation avec `geom_col` et `facet_wrap` ggplot(tab_piv, aes(x = libelle_groupe_culture,y = value)) +geom_col(fill ="lightblue", colour ="black", position ="dodge") +labs(title ="Surface par culture", x="Culture", y ="%", fill ="Secteur") +geom_text(aes(label = value), hjust =-0.3, size =8/.pt, colour ="black") +theme_classic() +coord_flip() +facet_wrap(~secteur, nrow=3, scales='free')```:::## Création d'un dashboard de visualisationOn souhaiterait intégrer les analyses faites ci-dessus à un tableau de bord qui offrirait des visualisations interactives. Une solution est d'utiliser `shiny`, une librairie qui permet la création de telles applications avec `R` notamment.Une application `shiny` simple a deux composants principaux :- Un objet d'interface utilisateur (UI) qui contrôle la disposition et l'apparence du tableau de bord ;- Une fonction serveur qui contient les instructions nécessaires au fonctionnement de l'application.Lorsque ces deux composants sont définis, l'application est lancée à l'aide d'une simple instruction `runApp`. Ici, on aimerait bien aller vers une première version de tableau de bord, qui afficherait au départ une carte interactive et laisserait aussi à l'utilisateur de choisir un rayon. Une fois ce rayon choisi et un point sélectionné grâce à un clic sur la carte, on souhaiterait (comme affiché ci-dessous) :- Afficher sur la carte les parcelles se situant à une distance inférieure au rayon choisi du point spécifié ;- Afficher des statistiques sur les parcelles en questions. Ainsi notre UI sera un objet [`fluidPage`](https://shiny.posit.co/r/reference/shiny/1.3.1/fluidpage) composé de trois éléments, et dans la fonction de serveur il y aura deux évènements distincts à [observer](https://shiny.posit.co/r/reference/shiny/1.0.1/observeevent). Pour développer l'application, créez un répertoire `my_app` et à l'intérieur 3 fichiers, un fichier `ui.R`, un fichier `server.R` et un fichier `utils.R` (dans lequel se trouveront les fonctions utilitaires).## UILe fichier UI aura la forme suivante :```{r}#| label: ui-q#| eval: false#| code-fold: falselibrary(shiny)# Define UIui <-fluidPage( ...)```A l'intérieur de la fonction `fluidPage` et grâce à la documentation de `Shiny`, ajoutez les 3 éléments qui constitueront le tableau de bord : un output de type "carte `leaflet`", un champ d'input numérique, et un output de type "table".::: {.callout-note icon=false}## Solution```{r}#| label: ui-a#| eval: falselibrary(shiny)# Define UIui <-fluidPage(leafletOutput("map", height =800),numericInput("buffer_radius", "Rayon (en km) :", value =5),tableOutput("table"))```:::## ServeurLe code de la fonction serveur figure ci-dessous dans son intégralité, et l'objectif de cette partie va être d'implémenter dans `utils.R` les différentes fonctions appelées dans le code :- `connect_to_db` : une fonction qui renvoie une connexion à la base de données ;- `query_db` : une fonction qui interroge la base de données pour récupérer les parcelles se situant à l'intérieur d'un cercle de rayon `radius` d'un point défini par une latitude et une longitude ;- `plot_surroundings` : une fonction qui prend en entrée la sortie de la fonction `query_db` et une carte `leaflet` et qui ajoute une couche comportant les polygones des parcelles concernées ;- `compute_stats` : une fonction qui prend en entrée la sortie de la fonction `query_db` et renvoie des statistiques sur les parcelles concernées.```{r}#| label: server-q#| eval: false#| code-fold: falselibrary(shiny)library(leaflet)# Définition du serveurserver <-function(input, output) {# Rendre la carte output$map <-renderLeaflet({leaflet() %>%addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%setView(lng =-1.4932577154775046, lat =46.46946179131805, zoom =12) })# Connexion à la base de données cnx <-connect_to_db()# Initialisation d'une "reactive value" pour le point sélectionné selectedPoint <-reactiveValues(lat =NULL, lng =NULL)# Gestion de l'évènement "clic"observeEvent(input$map_click, { clickData <- input$map_clickif (!is.null(clickData)) {# Stockage des coordonnées du point selectedPoint$lat <- clickData$lat selectedPoint$lng <- clickData$lng buffer_radius <- input$buffer_radius sf <-query_db(cnx, selectedPoint$lat, selectedPoint$lng, buffer_radius)# Mise à jour de la carteleafletProxy("map") %>%clearMarkers() %>%clearShapes() %>%addMarkers(lng = selectedPoint$lng, lat = selectedPoint$lat) %>%plot_surroundings(sf)# Calculs sur les données de parcelles df <-compute_stats(sf)# Update de la table affichée output$table <-renderTable({ df }) } })observeEvent(input$buffer_radius, {# Vérification qu'un point a été sélectionnéif (!is.null(selectedPoint$lat) &&!is.null(selectedPoint$lng)) {# Requête avec le nouveau rayon buffer_radius <- input$buffer_radius sf <-query_db(cnx, selectedPoint$lat, selectedPoint$lng, buffer_radius)# Mise à jour de la carteleafletProxy("map") %>%clearShapes() %>%plot_surroundings(sf)# Calculs sur les données de parcelles df <-compute_stats(sf)# Update de la table affichée output$table <-renderTable({ df }) } })}```Proposer une implémentation des fonctions `connect_to_db`, `query_db` et `compute_stats` décrites ci-dessus et utilisées dans la fonction "serveur".### `connect_to_db`La fonction `connect_to_db` renvoie une connexion à la base de données PostgreSQL, comme vu précédemment. Les identifiants (`user` et `password` sont respectivement stockés dans les variables d'environnement `USER_POSTGRESQL` et `PASS_POSTGRESQL`.)::: {.callout-note icon=false}## Solution```{r}#| label: connect-fun#| eval: false#' Connection au serveur PostgreSQL. Le mot de passe doit être stocké dans la#' variable d'environnement PASS_POSTGRESQL.#' #' @returns Connexion au serveur.connect_to_db <-function() {# Connection à PostgreSQL cnx <-dbConnect(Postgres(),user =Sys.getenv("USER_POSTGRESQL"),password =Sys.getenv("PASS_POSTGRESQL"),host ="postgresql-100400.projet-funathon",dbname ="defaultdb",port =5432,check_interrupts =TRUE)return(cnx)}```:::### `query_db`On souhaite implémenter la fonction suivante :```{r}#| label: query-fun-q#| eval: false#| code-fold: false#' Requête la table `parcelles` pour récupérer les parcelles qui se situent#' dans un certain rayon autour d'un point repéré par une latitude et #' une longitude.#' #' @param cnx Connexion à PostgreSQL.#' @param lat Latitude.#' @param lng Longitude.#' @param radius Rayon.#' @returns Objet `sf` avec les données des parcelles concernées.query_db <-function(cnx, lat, lng, radius) { ...}```On pourra notamment utiliser plusieurs fonctions PostGIS :- [`ST_MakePoint`](https://postgis.net/docs/ST_MakePoint.html) qui permet de créer une géométrie POINT;- [`ST_SetSRID`](https://postgis.net/docs/ST_SetSRID.html) qui permet de définir le système de coordonnées pour une géométrie;- [`ST_Buffer`](https://postgis.net/docs/ST_Buffer.html) qui calcule un POLYGON ou un MULTIPOLYGON qui représente tous les points dont la distance par rapport à une géométrie/géographie est inférieure ou égale à une distance donnée;- [`ST_Intersects`](https://postgis.net/docs/ST_Intersects.html) qui compare deux géométries et renvoie `true` si elles ont une intersection non-nulle;ainsi que des fonction de la librairie `sf`.::: {.callout-note icon=false}## Solution```{r}#| label: query-fun-a#| eval: false#' Requête la table `parcelles` pour récupérer les parcelles qui se situent#' dans un certain rayon autour d'un point repéré par une latitude et #' une longitude.#' #' @param cnx Connexion à PostgreSQL.#' @param lat Latitude.#' @param lng Longitude.#' @param radius Rayon.#' @returns Objet `sf` avec les données des parcelles concernées.query_db <-function(cnx, lat, lng, radius) {# Les données spatiales sur PostgreSQL sont stockées en Lambert 93.# Pour faire le join on veut donc projeter les coordonnées `lat`` et `lng` postgis_crs <-"EPSG:2154" coordinates <-data.frame(lng =c(lng), lat =c(lat)) %>%st_as_sf(coords =c("lng", "lat"), remove =TRUE) %>%st_set_crs("EPSG:4326") %>%st_transform(postgis_crs)# Requête PostgreSQL query <-sprintf("SELECT * FROM rpg.parcelles WHERE ST_Intersects(geom, ST_Buffer(ST_SetSRID(ST_MakePoint(%f, %f), 2154), %.0f));",st_coordinates(coordinates)[1],st_coordinates(coordinates)[2], radius*1000)# Récupération des résultats sf <-st_read( cnx,query = query )return(sf)}```:::### `compute_stats`On souhaite implémenter la fonction suivante :```{r}#| label: stats-fun-q#| eval: false#| code-fold: false#| #' Crée la table à afficher sur l'application grâce à des calculs sur les#' données requêtées depuis PostgreSQL.#' #' @param sf Données spatiales.#' @returns data.frame à afficher.compute_stats <-function(sf) { ...}```qui calcule sur une table spatiale renvoyée par la fonction `query_db` (ensemble des parcelles se situant autour d'un point donné) les statistiques suivantes par groupe de culture agregé :- Nombre de parcelles;- Pourcentage des parcelles parmi toutes les parcelles;- Surface des parcelles;- Surface des parcelles rapportée à la surface de toutes les parcelles;- Surface moyenne d'une parcelle.On récupère les groupes de cultures agrégés sur l'espace de stockage du SSP Cloud et avec un léger prétraitement, comme précédemment :```{r}#| label: lib-cult#| code-fold: false# Récupération des libellés des différentes cultureslib_cult <-s3read_using(FUN = read_csv2,object ="2023/sujet2/diffusion/ign/rpg/REF_CULTURES_GROUPES_CULTURES_2020.csv",col_types =cols(.default =col_character()),bucket ="projet-funathon",opts =list("region"="")) %>%clean_names()lib_group_cult <- lib_cult %>%select(code_groupe_culture, libelle_groupe_culture) %>%distinct(code_groupe_culture, .keep_all=T)lib_group_cult %>%kable()```::: {.callout-note icon=false}## Solution```{r}#| label: stats-fun-a#| eval: false#' Crée la table à afficher sur l'application grâce à des calculs sur les#' données requêtées depuis PostgreSQL.#' #' @param sf Données spatiales.#' @returns data.frame à afficher.compute_stats <-function(sf) { df <- sf %>%st_drop_geometry() %>%count(code_group, name ="parcelles_grp") %>%add_tally(parcelles_grp, name ="parcelles_tot") %>%mutate(pct_parcelles =round(100*parcelles_grp/parcelles_tot, 1)) %>%select(-parcelles_tot) %>%cbind(# Comptage des surfaces sf %>%st_drop_geometry() %>%count(code_group, wt = surf_parc, name ="surface_grp") %>%add_tally(surface_grp, name ="surface_tot") %>%mutate(surface_pct =round(100*surface_grp/surface_tot, 1)) %>%select(-surface_tot) %>%select(surface_grp, surface_pct) ) %>%left_join(lib_group_cult, by =c("code_group"="code_groupe_culture")) %>%select(code_group, libelle_groupe_culture, everything()) %>%arrange(desc(surface_grp)) %>%adorn_totals() %>%mutate(mean_surface =round(surface_grp/parcelles_grp, 1))return( df %>%select(-code_group) %>%setNames(c("Groupe de cultures","Nombre de parcelles","Pourcentage de parcelles","Surface (ha)","Surface (%)","Surface moyenne d'une parcelle (ha)")) )}```:::### `plot_surroundings`La fonction `plot_surroundings` prend en entrée une carte `leaflet` ainsi qu'un objet `sf` récupéré grâce à `query_db` et renvoie une carte où on a rajouté une couche avec les polygones correspondant aux parcelles. Elle est donné ci-dessous. Essayez de comprendre à quoi correspondent les différents arguments de la fonction `addPolygons`.```{r}#| label: plot-fun-a#| eval: false#| code-fold: false# Création d'une palette de couleurs associée au groupe de culturepal <-brewer.pal(12, "Paired")pal <-colorRampPalette(pal)(24)factpal <-colorFactor(pal, lib_group_cult$code_groupe_culture)#' Rajoute les données d'un objet `sf` sous forme de polygones à une#' carte `leaflet`.#' #' @param leaflet_proxy Carte.#' @param sf Données spatiales.#' @returns Carte enrichie.plot_surroundings <-function(leaflet_proxy, sf) {# Transformation de la projection (WGS 84) sf <- sf %>%st_transform(4326)# Ajout des libellés des cultures sf <- sf %>%left_join(lib_cult %>%select(-code_groupe_culture), by =c("code_cultu"="code_culture")) # Création des labels à afficher au passage de la souris sur la carte. labels <-sprintf("<strong>Identifiant de la parcelle : </strong>%s<br/> <strong>Groupe culture : </strong>%s<br/> <strong>Culture : </strong>%s<br/> <strong>Surface (ha) : </strong>%s<br/> <strong>Département : </strong>%s<br/> <strong>Commune : </strong>%s<br/>", sf$id_parcel, sf$libelle_groupe_culture, sf$libelle_culture, sf$surf_parc, sf$insee_dep, sf$nom_com) %>%lapply(htmltools::HTML)return( leaflet_proxy %>%addPolygons(data = sf,fillColor =~factpal(code_group),weight =2,opacity =1,color ="#ffd966",dashArray ="3",fillOpacity =0.5,highlight =highlightOptions(weight =5,color ="#A40000",dashArray ="",fillOpacity =0.0,bringToFront =TRUE),label = labels,labelOptions =labelOptions(style =list("font-weight"="normal", padding ="3px 8px"),textsize ="15px",direction ="auto",encoding="UTF-8")) )}```## Lancement de l'applicationL'application RShiny se lance localement grâce à la commande `shiny::runApp("my-app")`. Vérifiez que tout fonctionne bien ! Une application fonctionnelle se trouve dans le répertoire `app` si besoin.
Ainsi notre UI sera un objet