diff --git a/content/git/exogit.qmd b/content/git/exogit.qmd index d3b3a5a91..ae92623aa 100644 --- a/content/git/exogit.qmd +++ b/content/git/exogit.qmd @@ -145,11 +145,11 @@ les fonctionalités coeur de ces deux interfaces qui sont en fait quasi-identiqu Les deux premières étapes se font sur `Github`. - - -{{% box status="exercise" title="Exercice" icon="fab fa-github" %}} - -**Exercice 1 : Créer un compte Github** +::: {.cell .markdown} +```{=html} + +``` +::: ## Deuxième étape: créer un *token* (jeton) HTTPS diff --git a/content/manipulation/03_geopandas_TP.qmd b/content/manipulation/03_geopandas_TP.qmd index e272692c4..434e6badd 100644 --- a/content/manipulation/03_geopandas_TP.qmd +++ b/content/manipulation/03_geopandas_TP.qmd @@ -68,7 +68,7 @@ Illustration du principe des données spatiales (documentation de `sf`, l'équiv [^noteQGIS]: D'ailleurs, le logiciel de cartographie spécialisé QGIS, s'appuie sur `Python` pour les manipulations de données nécessaires avant de réaliser une carte. -Ce chapitre illustre à partir d’exemples pratiques certains principes centraux de l’analyse de données: +Ce chapitre illustre à partir d’exemples pratiques certains principes centraux de l’analyse de données : - Manipulations sur les attributs des jeux de données ; - Manipulations géométriques ; @@ -104,25 +104,38 @@ cellule de _notebook_ `Jupyter`, le code suivant est à exécuter: !pip install topojson ``` -Les instructions d'installation du package `cartiflette` -sont quant à elles détaillées dans le chapitre -précédent. - Après installations, les _packages_ à importer pour progresser dans ce chapitre sont les suivants: - ```{python} #| echo: true #| output: false import geopandas as gpd import contextily as ctx import matplotlib.pyplot as plt +``` + +Les instructions d'installation du package `cartiflette` +sont quant à elles détaillées dans le chapitre +précédent. + +```{python} +#| echo: false +#| output: false +!pip install requests py7zr geopandas openpyxl tqdm s3fs PyYAML xlrd +!pip install git+https://github.com/inseefrlab/cartiflette@80b8a5a28371feb6df31d55bcc2617948a5f9b1a +``` + +```{python} +#| echo: true +#| output: false from cartiflette.s3 import download_vectorfile_url_all ``` + + ## Lire et enrichir des données spatiales Dans cette partie, @@ -229,35 +242,33 @@ communes_borders.head() # Il y a une colonne geometry qui contient les informations nécessaires pour connaître les contours communaux ``` -A la question 3, le CRS à obtenir est: - ```{python} +#| output: false + # 3) Afficher le crs communes_borders.crs # Les données sont en WGS84, on les reprojette en lambert 93 ``` -La carte du 92 est la suivante, si on ne conserve que les limites des communes : - ```{python} +#| echo: false + # 4) afficher les communes du département 92 ax = communes_borders[communes_borders['INSEE_DEP'] == "92"].boundary.plot() ax.set_axis_off() ``` -Quant à Paris, à l'issue de la question 5, la même carte -aura l'aspect suivant: - ```{python} +#| echo: false + # 5) Représenter la carte de Paris. Quel est le problème ? ax = communes_borders[communes_borders['INSEE_DEP'] == "75"].boundary.plot() ax.set_axis_off() +#On remarque ainsi facilement le problème pour Paris +#intra-muros: il manque les limites des arrondissements. +#Cela appauvrit grandement la carte de Paris. ``` -On remarque ainsi facilement le problème pour Paris -intra-muros: il manque les limites des arrondissements. -Cela appauvrit grandement la carte de Paris. - ```{python} #| output: false @@ -271,17 +282,12 @@ petite_couronne = download_vectorfile_url_all( filter_by="DEPARTEMENT", source="EXPRESS-COG-CARTO-TERRITOIRE", year=2022) + +petite_couronne.crs +petite_couronne = petite_couronne.to_crs(2154) +petite_couronne.crs ``` -La carte de Paris intra-muros est, à l'issue de la -question 6, - -```{python} -ax = petite_couronne.boundary.plot(edgecolor = "black") -ax.set_axis_off() -``` - - ## Le système de projection Un concept central dans les logiciels de SIG est la notion de @@ -310,7 +316,7 @@ width_projected_map = screen.width/2 ::: {.cell .markdown} ```{=html} @@ -246,27 +450,31 @@ la récupération de fonds de carte sur un ensemble de département. Cela évite la récupération d'un fond de carte très volumineux (plus de 500Mo) pour une analyse restreinte (quelques départements). Un autre avantage de `cartiflette` est de faciliter la récupération de fonds -de carte consolidés comme celui dont on a besoin ici: arrondissements +de carte consolidés comme celui dont on a besoin ici : arrondissements dans Paris, communes ailleurs. Comme cela est expliqué dans un encadré à part, il s'agirait d'une opération pénible à mettre en oeuvre sans `cartiflette`. Les contours de cet espace peuvent être récupérés de la manière suivante : ```{python} +#| output: false import cartiflette.s3 as s3 shp_communes = s3.download_vectorfile_url_all( crs = 4326, values = ["75", "92", "93", "94"], - borders="COMMUNE", + borders="COMMUNE_ARRONDISSEMENT", vectorfile_format="topojson", filter_by="DEPARTEMENT", source="EXPRESS-COG-CARTO-TERRITOIRE", year=2022) +``` -shp_communes.head() +```{python} +shp_communes.head(3) ``` + On reconnaît la structure d'un `DataFrame` `Pandas`. A cette structure s'ajoute une colonne `geometry` qui enregistre la position des limites des polygones de chaque observation. @@ -308,6 +516,7 @@ carte en plusieurs phases. En premier lieu, il est nécessaire de récupérer le niveau des communes. ```{python} +#| output: false import cartiflette.s3 as s3 shp_communes = s3.download_vectorfile_url_all( @@ -318,8 +527,10 @@ shp_communes = s3.download_vectorfile_url_all( filter_by="DEPARTEMENT", source="EXPRESS-COG-CARTO-TERRITOIRE", year=2022) +``` -shp_communes.head() +```{python} +shp_communes.head(4) ``` ```{python} @@ -347,16 +558,15 @@ grâce à `cartiflette`: arrondissements = s3.download_vectorfile_url_all( crs = 4326, values = ["75"], - borders="COMMUNE", + borders="COMMUNE_ARRONDISSEMENT", vectorfile_format="topojson", filter_by="DEPARTEMENT", source="EXPRESS-COG-CARTO-TERRITOIRE", year=2022) - ``` ```{python} -ax = arrondissements.plot(alpha = 0.8, edgecolor = "k") +ax = arrondissements.boundary.plot(alpha = 0.8, edgecolor = "k") ax.set_axis_off() ``` @@ -378,6 +588,11 @@ shp_communes = pd.concat( ]) ``` +```{python} +ax = shp_communes.boundary.plot(alpha = 0.8, edgecolor = "k") +ax.set_axis_off() +``` + Cette approche fonctionne mais elle nécessite un certain nombre de gestes, qui sont autant de risques d'erreurs. Il est donc recommandé de privilégier le niveau `COMMUNE_ARRONDISSEMENT` @@ -390,7 +605,7 @@ qui fait exactement ceci mais de manière fiable. ## Opérations sur les attributs et les géométries -### Import des données velib +### Import des données vélib Souvent, le découpage communal ne sert qu'en fond de cartes, pour donner des repères. En complément de celui-ci, on peut désirer exploiter @@ -430,7 +645,7 @@ fig,ax = plt.subplots(figsize=(10, 10)) stations.sample(200).to_crs(3857).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2) shp_communes.to_crs(3857).plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none", color = None) -ctx.add_basemap(ax, source = ctx.providers.Stamen.Watercolor) +#ctx.add_basemap(ax, source = ctx.providers.Stamen.Watercolor) ax.set_axis_off() ``` @@ -444,25 +659,22 @@ fig.savefig("featured_geopandas_tutorial.png") Découvrez ci-dessous par étape les différentes lignes de commandes permettant d'afficher cette carte complète, étape par étape : -:one: -Afficher le nuage de points de 200 stations vélibs prises au hasard +1️⃣ Afficher le nuage de points de 200 stations vélibs prises au hasard ```{python} #| output: hide -fig,ax = plt.subplots(figsize=(10, 10)) +fig, ax = plt.subplots(figsize=(10, 10)) stations.sample(200).to_crs(3857).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2) ``` - -:two: -Ajouter à cette couche, en-dessous, les contours des communes +2️⃣ Ajouter à cette couche, en-dessous, les contours des communes ```{python} #| output: false #| echo: true -fig,ax = plt.subplots(figsize=(10, 10)) +fig, ax = plt.subplots(figsize=(10, 10)) stations.sample(200).to_crs(3857).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2) shp_communes.to_crs(3857).plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none", color = None) @@ -473,10 +685,12 @@ shp_communes.to_crs(3857).plot(ax = ax, zorder=1, edgecolor = "black", facecolor ax.get_figure() ``` -:three: -Ajouter un fond de carte de type *open street map* grâce au package + + +3️⃣ Ajouter un fond de carte de type *open street map* grâce au package `contextily` + ```{python} #| output: false #| echo: true @@ -485,7 +699,7 @@ fig,ax = plt.subplots(figsize=(10, 10)) stations.sample(200).to_crs(3857).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2) shp_communes.to_crs(3857).plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none", color = None) -ctx.add_basemap(ax, source = ctx.providers.Stamen.Watercolor) +#ctx.add_basemap(ax, source = ctx.providers.Stamen.Watercolor) ``` ```{python} @@ -494,7 +708,7 @@ ax.get_figure() ``` -:four: +4️⃣ Il ne reste plus qu'à retirer l'axe des coordonnées, qui n'est pas très esthétique. @@ -506,7 +720,7 @@ fig,ax = plt.subplots(figsize=(10, 10)) stations.sample(200).to_crs(3857).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2) shp_communes.to_crs(3857).plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none", color = None) -ctx.add_basemap(ax, source = ctx.providers.Stamen.Watercolor) +#ctx.add_basemap(ax, source = ctx.providers.Stamen.Watercolor) ax.set_axis_off() ax ``` @@ -526,51 +740,43 @@ sur un objet `GeoPandas`. Pour manipuler les données, et non la géométrie, on parlera d'opérations sur les attributs. Par exemple, si on désire -connaître quelques statistiques sur la taille des stations: +connaître quelques statistiques sur la taille des stations, l'approche +est identique à si on avait un objet `Pandas` classique : ```{python} stations.describe() ``` -Pour connaître les plus grands départements de France métropolitaine, +Pour classer les départements de la petite couronne, du plus grand au plus petit, procédons en deux étapes: -1. Récupérons le contour des communes de France métropolitaine dans son ensemble +1. Récupérons le contour des communes grâce à `cartiflette`. Notons qu'on pourrait récupérer directement les contours départementaux mais pour l'exercice, nous allons le créer nous-mêmes comme agrégation des contours communaux -(voir [ce notebook `observable`](https://observablehq.com/@linogaliana/cartiflette-demo) pour la méthode plus +(voir plus bas ainsi que [ce notebook `Observable`](https://observablehq.com/@linogaliana/cartiflette-demo) pour la méthode plus légère qui utilise pleinement les fonctionnalités de `cartiflette`). -2. Calculons la surface (méthode `area` sur un objet `GeoPandas.GeoDataFrame` ramenée en km², attention néamoins au système de projection comme cela est expliqué plus bas) -```{python} -#| output: false -from cartiflette.download import get_vectorfile_ign -france = get_vectorfile_ign( - borders = "COMMUNE", - field = "metropole", - source = "EXPRESS-COG-CARTO-TERRITOIRE", - provider="IGN" -) -``` +2. Calculons la surface totale de ce territoire (méthode `area` sur un objet `GeoPandas.GeoDataFrame` ramenée en km², attention néamoins au système de projection comme cela est expliqué plus bas) + ```{python} -france['surface'] = france.area.div(10**6) +shp_communes['surface'] = shp_communes.area.div(10**6) ``` Les plus grands départements s'obtiennent par une agrégation des surfaces communales : ```{python} -france.groupby('INSEE_DEP').sum(numeric_only = True).sort_values('surface', ascending = False) +shp_communes.groupby('INSEE_DEP').sum(numeric_only = True).sort_values('surface', ascending = False) ``` Si on veut directement les plus -grandes communes de France métropolitaine : +grandes communes de la petite couronne parisienne : ```{python} -france.sort_values('surface', ascending = False).head(10) +shp_communes.sort_values('surface', ascending = False).head(10) ``` Lors des étapes d'agrégation, `groupby` ne conserve pas les géométries. Autrement @@ -583,22 +789,38 @@ attribus avec la méthode `dissolve`: ```{python} fig,ax = plt.subplots(figsize=(10, 10)) -france.dissolve(by='INSEE_DEP', aggfunc='sum').plot(ax = ax, column = "surface") +shp_communes.dissolve(by='INSEE_DEP', aggfunc='sum').plot(ax = ax, column = "surface") ax.set_axis_off() ax ``` -Pour produire cette carte, il serait néanmoins plus simple de directement +Pour produire l'équivalent de cette carte à un niveau France entière, il est néanmoins plus simple de directement récupérer les fonds officiels des départements plutôt que d'agréger les contours des communes: ```{python} -dep = get_vectorfile_ign( - borders = "DEPARTEMENT", - field = "metropole", - source = "EXPRESS-COG-CARTO-TERRITOIRE", - provider="IGN") -dep["area"] = dep.area +#| output: false +dep = s3.download_vectorfile_url_all( + values = "metropole", + crs = 4326, + borders = "DEPARTEMENT", + vectorfile_format="topojson", + filter_by="FRANCE_ENTIERE", + source="EXPRESS-COG-CARTO-TERRITOIRE", + year=2022) + +dep["area"] = dep.to_crs(2154).area +``` + +Avant de calculer les surfaces des départements, pour éviter les déformations liées au +système `Mercator`, nous faisons une reprojection des données à la volée. Plus de détails +par la suite. + +```{python} +dep.sort_values('area', ascending = False).head(3) +``` + +```{python} ax = dep.plot(column = "area") ax.set_axis_off() ``` @@ -627,7 +849,7 @@ de carte `OpenStreetMap`, `Stamen`, `Google Maps`, etc. Mais ce n'est pas un format sur lequel on désire faire des calculs car les distances sont faussées sans utiliser de projection. D'ailleurs, `geopandas` refusera certaines opérations -sur des données dont le crs est `4326`. On reprojete ainsi les données +sur des données dont le crs est `4326`. On reprojette ainsi les données dans la projection officielle pour la métropole, le Lambert 93 (epsg: `2154`). @@ -681,7 +903,7 @@ communes_31['geometry'] = communes_31['geometry'].centroid ax = communes_31.plot(figsize = (10,10), color = 'red', alpha = 0.4, zorder=2) dep_31.to_crs(3857).plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none", color = None) -ctx.add_basemap(ax, source = ctx.providers.Stamen.Toner) +#ctx.add_basemap(ax, source = ctx.providers.Stamen.Toner) ax.set_axis_off() ax ``` @@ -699,37 +921,6 @@ stations = stations.to_crs(2154) ``` -Le système de projection est fondamental pour que la dimension -spatiale soit bien interprétée par `Python`. Un mauvais système de représentation -fausse l'appréciation visuelle mais peut aussi entraîner des erreurs dans -les calculs sur la dimension spatiale. -Ce [post](https://www.earthdatascience.org/courses/earth-analytics/spatial-data-r/geographic-vs-projected-coordinate-reference-systems-UTM/) propose de riches éléments sur le -sujet, notamment l'image suivante qui montre bien le principe d'une projection : - -![Les différents types de projection](https://www.earthdatascience.org/images/courses/earth-analytics/spatial-data/spatial-projection-transformations-crs.png) - -La Terre peut ainsi être représentée de multiples manières, ce qui n'est pas neutre sur la manière de se représenter -certains continents. -L'Afrique apparaît beaucoup moins vaste qu'elle ne l'est en réalité sur les cartes utilisant -cette projection. -L'une des déformations les mieux connue est celle provoquée par la projection Mercator. -Le Groënland paraît avoir -la même surface que l'Amérique du Sud. Pourtant, cette dernière est 8 fois -plus grande. - -Il existe en fait de nombreuses représentations possibles du monde, plus ou moins -alambiquées. Les projections sont très nombreuses et certaines peuvent avoir une [forme suprenante](https://imgs.xkcd.com/comics/map_projections.png). -Par exemple, -la [projection de Spillhaus](https://storymaps.arcgis.com/stories/756bcae18d304a1eac140f19f4d5cb3d) -propose de centrer la vue sur les océans et non une terre. C'est pour -cette raison qu'on parle parfois de monde tel que vu par les poissons -à son propos. - - -![](truesize.png) -*Exemple de reprojection de pays depuis le site [thetruesize.com](https://thetruesize.com/)* - - Concernant la gestion des projections avec `GeoPandas`, la [documentation officielle](https://geopandas.org/projections.html) est très bien faite. Elle fournit notamment l'avertissement suivant qu'il est @@ -739,8 +930,8 @@ bon d'avoir en tête: > > From time to time, however, you may get data that does not include a projection. In this situation, you have to set the CRS so geopandas knows how to interpret the coordinates. -![](https://imgs.xkcd.com/comics/bad_map_projection_south_america.png) -*Image empruntée à XKCD qu'on peut également trouver sur * +![*Image empruntée à XKCD qu'on peut également trouver sur *](https://imgs.xkcd.com/comics/bad_map_projection_south_america.png) + Pour déterminer le système de projection d'une base de données, on peut vérifier l'attribut `crs`: @@ -748,7 +939,7 @@ Pour déterminer le système de projection d'une base de données, on peut véri communes.crs ``` -Les deux principales méthodes pour définir le système de projection utilisé sont: +Les deux principales méthodes pour définir le système de projection utilisé sont : * **`df.set_crs`**: cette commande sert à préciser quel est le système de projection utilisé, c'est-à-dire comment les coordonnées *(x,y)* sont reliées à la surface terrestre. **Cette commande ne doit pas être utilisée pour transformer le système de coordonnées, seulement pour le définir**. * **`df.to_crs`**: **cette commande sert à projeter les points d'une géométrie dans une autre, c'est-à-dire à recalculer les coordonnées selon un autre système de projection.** @@ -764,10 +955,15 @@ communes = communes.set_crs(2154) Alors que la reprojection (projection Albers : `5070`) s'obtient de la manière suivante : ```{python} -shp_region = get_vectorfile_ign( - borders = "REGION", field = "metropole", - source = "EXPRESS-COG-CARTO-TERRITOIRE", provider="IGN") - +shp_region = s3.download_vectorfile_url_all( + values = "metropole", + crs = 4326, + borders = "REGION", + vectorfile_format="topojson", + filter_by="FRANCE_ENTIERE", + source="EXPRESS-COG-CARTO-TERRITOIRE", + year=2022) + fig,ax = plt.subplots(figsize=(10, 10)) shp_region.to_crs(5070).plot(ax = ax) ax @@ -863,120 +1059,3 @@ Les jointures spatiales peuvent être très gourmandes en ressources (car il peu ``` ::: -## Annexe - -### Récupération des données depuis data.gouv - -Avec cette méthode, les données des limites administratives demandent donc un peu de travail pour être -importées car elles sont zippées (mais c'est un bon exercice !). - -Le code suivant, dont les -détails apparaîtront plus clairs après la lecture de la partie -*[webscraping](webscraping)* permet de : - -1. Télécharger les données avec `requests` dans un dossier temporaire -2. Les dézipper avec le module `zipfile` - -La fonction suivante automatise un peu le processus : - -```{python} -#| eval: false -import requests -import tempfile -import zipfile - -url = 'https://www.data.gouv.fr/fr/datasets/r/0e117c06-248f-45e5-8945-0e79d9136165' -temporary_location = tempfile.gettempdir() - -def download_unzip(url, dirname = tempfile.gettempdir(), destname = "borders"): - myfile = requests.get(url) - open("{}/{}.zip".format(dirname, destname), 'wb').write(myfile.content) - with zipfile.ZipFile("{}/{}.zip".format(dirname, destname), 'r') as zip_ref: - zip_ref.extractall(dirname + '/' + destname) -``` - -```{python} -#| eval: false -#| output: false - -download_unzip(url) -``` - -```{python} -#| eval: false -#| output: false -shp_communes = gpd.read_file(temporary_location + "/borders/communes-20220101.shp") -``` - -Ici, les données ne sont pas projetées puisqu'elles sont dans le -système WSG84 (epsg: 4326) ce qui permet de facilement ajouter -un fonds de carte `Openstreetmap` ou `Stamen` pour rendre une représentation -graphique plus esthétique. -En toute rigueur, pour faire une carte statique d'un pays en particulier, -il faudrait reprojeter les données dans un système de projection adapté à la zone géographique étudiée -(par exemple le Lambert-93 pour la France métropolitaine). - -On peut ainsi représenter Paris pour se donner une idée de la nature -du shapefile utilisé : - -```{python} -#| eval: false -paris = shp_communes.loc[shp_communes['insee'].str.startswith("75")] - -fig,ax = plt.subplots(figsize=(10, 10)) -paris.plot(ax = ax, alpha=0.5, edgecolor='blue') -ctx.add_basemap(ax, crs = paris.crs.to_string()) -ax.set_axis_off() -ax -``` - -On voit ainsi que les données pour Paris ne comportent pas d'arrondissement, -ce qui est limitant pour une analyse focalisée sur Paris. On va donc les -récupérer sur le site d'open data de la ville de Paris et les substituer -à Paris : - -```{python} -#| eval: false -#| echo: true - -arrondissements = gpd.read_file("https://opendata.paris.fr/explore/dataset/arrondissements/download/?format=geojson&timezone=Europe/Berlin&lang=fr") -arrondissements = arrondissements.rename(columns = {"c_arinsee": "insee"}) -arrondissements['insee'] = arrondissements['insee'].astype(str) -shp_communes = shp_communes[~shp_communes.insee.str.startswith("75")].append(arrondissements) -``` - -Pour produire la carte, il faudrait faire: - -```{python} -#| eval: false -paris = shp_communes.loc[shp_communes.insee.str.startswith("75")] - -fig,ax = plt.subplots(figsize=(10, 10)) - -paris.plot(ax = ax, alpha=0.5, edgecolor='k') -ctx.add_basemap(ax, crs = paris.crs.to_string()) -ax.set_axis_off() -ax -``` - - -### Récupération des données depuis le package `pynsee` - - -Pour connaître les contraintes d'installation du package `pynsee`, -se référer à la partie de cours dédiée à `Pandas`. - -```{python} -#| output: false -#| eval: false - -#le téléchargement des données prend plusieurs minutes -from pynsee.geodata import get_geodata -shp_communes = gpd.GeoDataFrame( - get_geodata('ADMINEXPRESS-COG-CARTO.LATEST:commune') -) -shp_communes = shp_communes.rename({"insee_com": 'insee'}, axis = 'columns') -#shp_communes = shp_communes.set_crs(3857) -``` - - diff --git a/content/manipulation/_play_with_projection.qmd b/content/manipulation/_play_with_projection.qmd index ac341693f..08ce19ead 100644 --- a/content/manipulation/_play_with_projection.qmd +++ b/content/manipulation/_play_with_projection.qmd @@ -1,4 +1,5 @@ ```{ojs} +//| echo: false container_projection = html`
@@ -16,6 +17,7 @@ container_projection = html`
```{ojs} +//| echo: false viewof projection = projectionInput({ name: "", value: "Mercator" @@ -23,11 +25,13 @@ viewof projection = projectionInput({ ``` ```{ojs} +//| echo: false import {projectionInput} from "@fil/d3-projections" import {map} from "@linogaliana/base-map" ``` ```{ojs} +//| echo: false projectedMap = map(projection, { //svg: true,