Finalize geopandas section (#48)

Encore le problème des figures (#29)

Author: linogaliana

config.toml

@@ -6,6 +6,8 @@ theme="hugo-theme-techdoc"
 enableGitInfo = true
 enableEmoji = true
 
+ignoreFiles = ["\\.Rmd$", "\\.Rmarkdown$", "_files$", "_cache$"]
+
 disableToc = false
 
 [markup.goldmark.renderer]

content/manipulation/03_geopandas.Rmd

@@ -25,35 +25,52 @@ slug: geopandas
 library(knitr)  
 library(reticulate)  
 knitr::knit_engines$set(python = reticulate::eng_python)
+# knitr::opts_knit$set(#fig.path = '../../static/Rmdoutput/manipulation',
+#                      root.dir = '../../static/Rmdoutput/manipulation'#,
+#                      #base.dir = '../../static/Rmdoutput/manipulation'#,
+#                      #output.dir = '../../static/Rmdoutput/manipulation'
+#                      )
 ```
 
 ```{python, include = FALSE}
 import os
 os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = 'C:/Users/W3CRK9/AppData/Local/r-miniconda/envs/r-reticulate/Library/plugins/platforms'
+os.environ["PROJ_LIB"] = r'C:\Users\W3CRK9\AppData\Local\r-miniconda\pkgs\proj4-4.9.3-hfa6e2cd_9\Library\share'
+#os.environ['GDAL_DATA'] = "C:/Users/W3CRK9/AppData/Local/gdal-3-1-2/bin/gdal-data"
 ```
 
 Dans ce tutoriel, nous allons utiliser:
 
 <!----
 * [Réseau de pistes cyclables de la ville de Paris](https://opendata.paris.fr/explore/dataset/reseau-cyclable/table/?disjunctive.typologie_simple&disjunctive.bidirectionnel&disjunctive.statut&disjunctive.sens_velo&disjunctive.voie&disjunctive.arrdt&disjunctive.bois&disjunctive.position&disjunctive.circulation&disjunctive.piste&disjunctive.couloir_bus&disjunctive.type_continuite&disjunctive.reseau&basemap=jawg.streets&location=12,48.85943,2.3493)
 ------>
-* [Carte des limites administratives françaises]()
+
+* [Localisations des stations velib](https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
+* [Carte des limites administratives françaises](https://www.data.gouv.fr/fr/datasets/r/07b7c9a2-d1e2-4da6-9f20-01a7b72d4b12)
+* [Arrondissements parisiens](https://opendata.paris.fr/explore/dataset/arrondissements/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
 
 La représentation des données, notamment la cartographie, est présentée plus
 amplement dans la partie [visualiser](visualiser). 
 
-Ce tutoriel s'inspire beaucoup d'un autre tutoriel que j'ai fait pour R **lien utilitR**. Il peut servir de pendant à celui-ci pour l'utilisateur de `R`. 
+Ce tutoriel s'inspire beaucoup d'un autre tutoriel que j'ai fait pour
+R disponible [ici](https://linogaliana.gitlab.io/documentationR/spatdata.html).
+Il peut servir de pendant à celui-ci pour l'utilisateur de `R`. 
 
 
 ```{python}
 import geopandas as gpd
-#import contextily as ctx
+import contextily as ctx
 ```
 
 
 ## Données spatiales: quelle différence avec des données traditionnelles ?
 
-**Le terme "données spatiales" désigne les données qui portent sur les caractéristiques géographiques des objets (localisation, contours, liens)**. Les caractéristiques géographiques des objets sont décrites à l'aide d'un **système de coordonnées** qui permetent une représentation dans un espace euclidien ($(x,y)$). Le passage de l'espace réel (la Terre, qui est une sphère) à l'espace plan se fait grâce à un **système de projection**. Voici quelques exemples de données spatiales:
+**Le terme "données spatiales" désigne les données
+qui portent sur les caractéristiques géographiques des objets (localisation, contours, liens)**. Les caractéristiques géographiques des objets sont décrites à l'aide d'un **système de coordonnées**
+qui permettent une représentation dans un espace euclidien ($(x,y)$).
+Le passage de l'espace réel (la Terre, qui est une sphère) à l'espace plan
+se fait grâce à un **système de projection**. Voici quelques exemples
+de données spatiales:
 
 * Une table décrivant des bâtiments, avec les coordonnées géographiques de chaque bâtiment;
 * Le découpage communal du territoire, avec le contour du territoire de chaque commune;
@@ -243,7 +260,7 @@ communes['dep'] = communes.insee.str[:2]
 ax = stations.sample(200).plot(figsize = (10,10), color = 'red', alpha = 0.4, zorder=2)
 communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
                                                            color = None)
-# ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
+#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
 ```
 
 ### Opérations sur les attributs
@@ -324,7 +341,7 @@ departement['geometry'] = departement['geometry'].centroid
 ax = departement.plot(figsize = (10,10), color = 'red', alpha = 0.4, zorder=2)
 communes[communes['dep'] == "31"].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
                                                            color = None)
-#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = #ctx.providers.Stamen.Toner)
+#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Toner)
 ```
 
 
@@ -334,16 +351,18 @@ Précédemment, nous avons appliqué une méthode `to_crs` pour reprojeter
 les données dans un système de projection différent de celui du fichier
 d'origine:
 
-```{python, eval = FALSE}
+```{python, eval = TRUE}
 communes = communes.to_crs(2154)
 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
+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:
+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)
 
@@ -377,15 +396,15 @@ communes = communes.set_crs(2154)
 
 Alors que la reprojection (projection Albers: 5070) s'obtient de la manière suivante:
 
-```{python, eval = FALSE}
+```{python, eval = TRUE}
 communes[communes.dep != "97"].dissolve(by='dep', aggfunc='sum').to_crs(5070).plot()
 ```
 
 On le voit, cela modifie totalement la représentation de l'objet dans l'espace.
 Clairement, cette projection n'est pas adaptée aux longitudes et lattitudes françaises. C'est normal, il s'agit d'une projection adaptée au continent 
 nord-américain, et encore, pas dans son ensemble :
 
-```{python, eval = FALSE}
+```{python, eval = TRUE}
 world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
 world[world.continent == "North America"].to_crs(5070).plot(alpha = 0.2, edgecolor = "k")
 ```
@@ -403,16 +422,16 @@ La seule différence avec `pandas` est dans la dimension géographique.
 Si on désire conserver la dimension géographique, il faut faire 
 attention à faire:
 
-```python
+~~~~python
 geopandas_object.merge(pandas_object)
-```
+~~~~
 
 Si on utilise deux objets géographiques mais ne désire conserver qu'une seule
 dimension géographique^[1], on fera
 
-```python
+~~~~python
 geopandas_object1.merge(geopandas_object2)
-```
+~~~~
 
 Seule la géométrie de l'objet de gauche
 sera conservée, même si on fait un *right join*. 

content/manipulation/03_geopandas.md

@@ -30,24 +30,34 @@ Dans ce tutoriel, nous allons utiliser:
 <!----
 * [Réseau de pistes cyclables de la ville de Paris](https://opendata.paris.fr/explore/dataset/reseau-cyclable/table/?disjunctive.typologie_simple&disjunctive.bidirectionnel&disjunctive.statut&disjunctive.sens_velo&disjunctive.voie&disjunctive.arrdt&disjunctive.bois&disjunctive.position&disjunctive.circulation&disjunctive.piste&disjunctive.couloir_bus&disjunctive.type_continuite&disjunctive.reseau&basemap=jawg.streets&location=12,48.85943,2.3493)
 ------>
-* [Carte des limites administratives françaises]()
+
+* [Localisations des stations velib](https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
+* [Carte des limites administratives françaises](https://www.data.gouv.fr/fr/datasets/r/07b7c9a2-d1e2-4da6-9f20-01a7b72d4b12)
+* [Arrondissements parisiens](https://opendata.paris.fr/explore/dataset/arrondissements/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
 
 La représentation des données, notamment la cartographie, est présentée plus
 amplement dans la partie [visualiser](visualiser). 
 
-Ce tutoriel s'inspire beaucoup d'un autre tutoriel que j'ai fait pour R **lien utilitR**. Il peut servir de pendant à celui-ci pour l'utilisateur de `R`. 
+Ce tutoriel s'inspire beaucoup d'un autre tutoriel que j'ai fait pour
+R disponible [ici](https://linogaliana.gitlab.io/documentationR/spatdata.html).
+Il peut servir de pendant à celui-ci pour l'utilisateur de `R`. 
 
 
 
 ```python
 import geopandas as gpd
-#import contextily as ctx
+import contextily as ctx
 ```
 
 
 ## Données spatiales: quelle différence avec des données traditionnelles ?
 
-**Le terme "données spatiales" désigne les données qui portent sur les caractéristiques géographiques des objets (localisation, contours, liens)**. Les caractéristiques géographiques des objets sont décrites à l'aide d'un **système de coordonnées** qui permetent une représentation dans un espace euclidien ($(x,y)$). Le passage de l'espace réel (la Terre, qui est une sphère) à l'espace plan se fait grâce à un **système de projection**. Voici quelques exemples de données spatiales:
+**Le terme "données spatiales" désigne les données
+qui portent sur les caractéristiques géographiques des objets (localisation, contours, liens)**. Les caractéristiques géographiques des objets sont décrites à l'aide d'un **système de coordonnées**
+qui permettent une représentation dans un espace euclidien ($(x,y)$).
+Le passage de l'espace réel (la Terre, qui est une sphère) à l'espace plan
+se fait grâce à un **système de projection**. Voici quelques exemples
+de données spatiales:
 
 * Une table décrivant des bâtiments, avec les coordonnées géographiques de chaque bâtiment;
 * Le découpage communal du territoire, avec le contour du territoire de chaque commune;
@@ -259,7 +269,7 @@ communes['dep'] = communes.insee.str[:2]
 ax = stations.sample(200).plot(figsize = (10,10), color = 'red', alpha = 0.4, zorder=2)
 communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
                                                            color = None)
-# ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
+#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
 ```
 
 ![](03_geopandas_files/figure-html/unnamed-chunk-5-1.png)<!-- -->
@@ -277,9 +287,9 @@ stations.describe()
 
 ```
 ##           capacity
-## count  1397.000000
-## mean     31.404438
-## std      12.149572
+## count  1398.000000
+## mean     31.403433
+## std      12.145281
 ## min       0.000000
 ## 25%      23.000000
 ## 50%      29.000000
@@ -327,11 +337,11 @@ communes[communes.dep != "97"].sort_values('surf_ha', ascending = False)
 ## 32467  49092         Chemillé-en-Anjou  ...   NaN   49
 ## 28877  49228           Noyant-Villages  ...   NaN   49
 ## ...      ...                       ...  ...   ...  ...
-## 15     75113                       NaN  ...  13.0   75
-## 16     75109                       NaN  ...   9.0   75
-## 17     75118                       NaN  ...  18.0   75
-## 18     75111                       NaN  ...  11.0   75
-## 19     75116                       NaN  ...  16.0   75
+## 15     75101                       NaN  ...   1.0   75
+## 16     75102                       NaN  ...   2.0   75
+## 17     75119                       NaN  ...  19.0   75
+## 18     75115                       NaN  ...  15.0   75
+## 19     75112                       NaN  ...  12.0   75
 ## 
 ## [34870 rows x 13 columns]
 ```
@@ -394,11 +404,17 @@ fera
 departement = communes[communes.dep == "31"].copy()
 departement['geometry'] = departement['geometry'].centroid
 
+```
+
+```
+## C:/Users/W3CRK9/AppData/Local/r-miniconda/envs/r-reticulate/python.exe:1: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
+```
 
+```python
 ax = departement.plot(figsize = (10,10), color = 'red', alpha = 0.4, zorder=2)
 communes[communes['dep'] == "31"].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
                                                            color = None)
-#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = #ctx.providers.Stamen.Toner)
+#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Toner)
 ```
 
 ![](03_geopandas_files/figure-html/unnamed-chunk-11-1.png)<!-- -->
@@ -417,10 +433,12 @@ 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
+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:
+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)
 
@@ -443,7 +461,20 @@ communes.crs
 ```
 
 ```
-## {'init': 'epsg:4326'}
+## <Projected CRS: EPSG:2154>
+## Name: RGF93 / Lambert-93
+## Axis Info [cartesian]:
+## - X[east]: Easting (metre)
+## - Y[north]: Northing (metre)
+## Area of Use:
+## - name: France
+## - bounds: (-9.86, 41.15, 10.38, 51.56)
+## Coordinate Operation:
+## - name: Lambert-93
+## - method: Lambert Conic Conformal (2SP)
+## Datum: Reseau Geodesique Francais 1993
+## - Ellipsoid: GRS 1980
+## - Prime Meridian: Greenwich
 ```
 
 Les deux principales méthodes pour définir le système de projection utilisé sont:
@@ -465,6 +496,8 @@ Alors que la reprojection (projection Albers: 5070) s'obtient de la manière sui
 communes[communes.dep != "97"].dissolve(by='dep', aggfunc='sum').to_crs(5070).plot()
 ```
 
+![](03_geopandas_files/figure-html/unnamed-chunk-15-1.png)<!-- -->
+
 On le voit, cela modifie totalement la représentation de l'objet dans l'espace.
 Clairement, cette projection n'est pas adaptée aux longitudes et lattitudes françaises. C'est normal, il s'agit d'une projection adaptée au continent 
 nord-américain, et encore, pas dans son ensemble :
@@ -475,6 +508,8 @@ world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
 world[world.continent == "North America"].to_crs(5070).plot(alpha = 0.2, edgecolor = "k")
 ```
 
+![](03_geopandas_files/figure-html/unnamed-chunk-16-1.png)<!-- -->
+
 
 ## Joindre des données
 
@@ -488,16 +523,16 @@ La seule différence avec `pandas` est dans la dimension géographique.
 Si on désire conserver la dimension géographique, il faut faire 
 attention à faire:
 
-```python
+~~~~python
 geopandas_object.merge(pandas_object)
-```
+~~~~
 
 Si on utilise deux objets géographiques mais ne désire conserver qu'une seule
 dimension géographique^[1], on fera
 
-```python
+~~~~python
 geopandas_object1.merge(geopandas_object2)
-```
+~~~~
 
 Seule la géométrie de l'objet de gauche
 sera conservée, même si on fait un *right join*. 

content/manipulation/03_geopandas_tp.Rmd

@@ -31,8 +31,13 @@ knitr::opts_chunk$set(eval = FALSE, include = FALSE, echo = FALSE)
 ```{python, include = FALSE}
 import os
 os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = 'C:/Users/W3CRK9/AppData/Local/r-miniconda/envs/r-reticulate/Library/plugins/platforms'
+os.environ["PROJ_LIB"] = r'C:\Users\W3CRK9\AppData\Local\r-miniconda\pkgs\proj4-4.9.3-hfa6e2cd_9\Library\share'
+#os.environ['GDAL_DATA'] = "C:/Users/W3CRK9/AppData/Local/gdal-3-1-2/bin/gdal-data"
 ```
 
+```{python, echo = TRUE}
+import geopandas as gpd
+```
 
 
 # Lire et enrichir des données spatiales
@@ -147,7 +152,6 @@ Cette [page](https://geopandas.org/mapping.html#maps-with-layers) peut vous aide
 ```{python}
 url = "https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr"
 stations = gpd.read_file(url)
-
 ```
 
 ```{python}
@@ -205,14 +209,33 @@ data_paris.merge(stations_agg, how = 'inner', suffixes = [None, '_agg']).plot(co
 ```{python}
 df = data_paris.merge(stations_agg, how = 'inner', suffixes = [None, '_agg'])
 cols = ['stationcode','capacity']
-df = df[[s + '_density' for s in cols]] = df[cols].div(df.to_crs(2158).area*10**(-6), axis = 0)
+df[[s + '_density' for s in cols]] = df[cols].div(df.to_crs(2158).area*10**(-6), axis = 0)
 df.plot(column = 'capacity_density', cmap = 'RdYlBu_r')
 df.plot(column = 'capacity_density', cmap = 'plasma_r')
 ```
 
 **Exercice 5 (optionnel): Relier distance au métro et capacité d'une station**
 
-1. Relier chaque station velib à la station de transport en commun la plus proche
+Une aide [ici](https://pysal.org/scipy2019-intermediate-gds/deterministic/gds1-relations.html#how-about-nearest-neighbor-joins)
+
+1. Relier chaque station velib à la station de transport en commun la plus proche. Vous pouvez
+prendre les localisations des stations [ici](https://data.iledefrance-mobilites.fr/explore/dataset/traces-du-reseau-ferre-idf/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
 2. Quelle ligne de transport est à proximité du plus de velib ?
 3. Calculer la distance de chaque station à la ligne de métro la plus proche. Faire un nuage de points reliant distance au métro et nombre de places en stations
 
+```{python}
+stations_transport = gpd.read_file("https://data.ratp.fr/explore/dataset/positions-geographiques-des-stations-du-reseau-ratp/download/?format=geojson&timezone=Europe/Berlin&lang=fr")
+```
+
+```{python, eval = FALSE}
+# un truc comme ça à revoir mais là fait des NaN :'(
+from shapely.ops import nearest_points
+temp = stations_transport.unary_union
+
+def find_nearest(geometry1, geometry2 = temp, df_point2 = stations_transport):
+    nearest = nearest_points(geometry1, geometry2)[1]
+    return df_point2[df_point2['geometry'] == nearest][:1]['nomptar']
+
+stations['Nearest'] = stations.apply(lambda row: find_nearest(row.geometry), axis=1)
+```
+