Manipuler des données géospatiales en Python¶
Les bibliothèques classiques permettant de manipuler des données géospatiales s'intègrent très bien aux IDE scientifiques (Spyder, Jupyter, etc.) et permettent de visualiser rapidement et interactivement des géométries ou un extrait d'une couche de données (attributs, géométries).
Les principales bibliothèques à connaitre sont les suivantes :
- shapely : bindings GEOS pour Python (équivalent R : rgeos)
- geopandas : manipulation de données vecteur sous la forme lignes/colonnes (avec pandas) (équivalent R : sf)
- rasterio : manipulation de données raster (équivalent R : raster)
- Rtree : création/interrogation d'un index spatial, bindings de libspatialindex (équivalent R : rtree de rgeos )
D'autres bibliothèques géospatiales sont également utilisables en Python :
Se familiariser avec les types géométriques de base : Shapely¶
La bibliothèque Shapely se présente comme des bindings de haut niveau autour de la bibliothèque GEOS.
In [3]:
from shapely.geometry import Point, LineString, Polygon
In [4]:
#### Création d'un point à partir de coordonnées :
pt1 = Point(1.0, 1.0)
pt2 = Point(3.1, 2.4)
# On peut utiliser 3 coordonnées :
pt_3d = Point(1.0, 1.0, 3.0)
# Accéder aux coordonnées d'un objet de type Point :
coords = pt1.coords[0]
print("Coordonnées sous forme d'un tuple avec pt.coords[0] : {}".format(coords))
# Accéder seulement à la coordonnée en 'x' ou en 'y' :
pt_x = pt1.x
pt_y = pt1.y
print("pt.x -> {}\npt.y -> {}\n".format(pt_x, pt_y))
# La représentation de ces geometry (lors de l'appel de print par exemple)
# utilise la notation WKT :
print(pt1)
In [5]:
# Dans Spyder comme dans Jupyter on pourra afficher une représentation graphique
# de la géométrie en inscrivant le nom de la variable qui la contient dans l'interpréteur :
pt1
Out[5]:
In [6]:
# Création d'une ligne à partir de plusieurs Point :
line = LineString([pt1, pt2])
# Création d'une ligne à partir de coordonnées :
line = LineString([(1.0, 1.0), (3.1, 2.4), (5.8, 4.5)])
# Accéder à la séquence de coordonnées de la ligne
# (sous forme d'une liste de tuples) :
list_coords = list(line.coords)
print("Sous forme d'une liste de tuples :\n{}\n".format(list_coords))
# Accéder seulement aux valeurs en 'x' :
list_xs = list(line.coords.xy[0])
print("Récupération des 'x' de notre séquence de points :\n{}\n".format(list_xs))
# Accéder seulement aux valeurs en 'y' :
list_ys = list(line.coords.xy[1])
print("Récupération des 'y' de notre séquence de points :\n{}\n".format(list_ys))
In [7]:
# Création d'un polygone à partir d'un seul anneau
# (l'anneau externe)
poly1 = Polygon([(1., 1.), (1.2, 1.), (1.2, 1.2), (1., 1.2), (1., 1.)])
# Création d'un polygone à partir de plusieurs anneaux
# (l'anneau externe et un anneau représentant un trou)
poly_with_hole = Polygon(
[(10., 10.), (12., 10.), (12., 12.), (10., 12.), (10., 10.)],
[[(11., 10.5), (10.5, 11.5), (11.5, 11.5), (11., 10.5)]])
print('Polygon:')
poly1
print('Polygon with hole:')
poly_with_hole
Out[7]:
Out[7]:
In [8]:
# Imports relatifs aux géométries
# multi-parties
from shapely.geometry import (
MultiPoint, MultiLineString,
MultiPolygon, GeometryCollection,
)
In [9]:
multipoint1 = MultiPoint([
Point(1.1, 1.2),
Point(8.0, 2.0),
Point(9.5, 9.5),
])
multipoint2 = MultiPoint([(1.1, 1.2), (8.0, 2.0), (9.5, 9.5)])
multiline = MultiLineString([
LineString([(12.1, 5.6), (4.5, 3.0)]),
LineString([(117.1, 12.6), (30.5, 46.0)]),
LineString([(0., 0.), (3., 3.)])
])
multipoly = MultiPolygon([poly1.buffer(2, 1), poly_with_hole])
collec = GeometryCollection([
Point(5., 5.),
LineString([(4.2, 4.3), (3., 3.), (5., 3.), (5., 0.,)]),
multipoly,
])
collec
Out[9]:
Les méthodes associées à ces objets :¶
In [10]:
# Des méthodes correspondants à des prédicats spatiaux :
poly1.contains(pt1)
multiline.intersects(poly1)
multipoint1.equals(multipoint2)
# Des propriétés :
poly1.boundary.is_closed
poly1.is_closed
poly1.area
poly1.intersects(poly_with_hole)
poly_with_hole.area
Out[10]:
Out[10]:
Out[10]:
Out[10]:
Out[10]:
Out[10]:
Out[10]:
Out[10]:
In [11]:
# Des méthodes qui permettent d'obtenir de nouvelles géométries :
pt1.buffer(10)
multiline.intersection(poly1)
# Des... propriétés qui permettent d'obtenir de nouvelles géométries :
multipoint1.convex_hull
poly_with_hole.boundary
Out[11]:
Out[11]:
Out[11]:
Out[11]:
Ce sont ces objets qu'on va manipuler dans la colonne geometry avec geopandas par la suite.
On retrouvera les attributs de ces objets sur cette colonne, avec pour signification d'appliquer l'opération à la géométrie de chaque entité.
geopandas¶
Ouverture d'une couche de données avec geopandas :
In [12]:
# Import nécessaires:
# Geopandas et les bibliothèques vues précédemment :
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
gdf = gpd.read_file('data/COURS_D_EAU.shp')
# L'objet est de type `GeoDataFrame` et
# possède les méthodes des `DataFrame` ainsi que d'autres :
gdf.head()
Out[12]:
In [13]:
# Attribut propre aux `GeoDataFrame` :
gdf.crs
Out[13]:
In [14]:
# Attributs dont le comportement a été modifié par rapport aux DataFrame :
gdf.plot()
Out[14]:
On dispose d'une colonne geometry. Lors de l'affichage de la table comme précédemment on a vu sa représentation au format WKT mais en interne il s'agit de géometries de la bibliothèque shapely.
Comme avec pandas, on peut sélectionner les colonnes avec la notation entre crochets et sélectionner des entités avec les méthodes .loc, .iloc ou en vérifiant une condition.
In [15]:
# Affichage de seulement quelques colonnes :
gdf[['TOPONYME', 'IMPORTANCE']].head()
Out[15]:
In [16]:
# Affichage de seulement quelques colonnes
# pour les cours d'eau dont l'importance est supérieure ou égale à 4:
gdf[gdf.IMPORTANCE >= 4][['TOPONYME', 'geometry']].head()
Lecture du message d'erreur et analyse du pbm ...¶
In [17]:
# Conversion explicite de type
gdf.IMPORTANCE = gdf.IMPORTANCE.astype(int)
# Retour à l'opération demandée
gdf[gdf.IMPORTANCE >= 4][['TOPONYME', 'geometry']].head()
Out[17]:
In [18]:
# Obtention d'une Serie contenant les longueurs de cours d'eau
lengths = gdf.geometry.length
lengths.sort_values(ascending=False).head() # <- tri effectué et ID conservés
Out[18]:
In [19]:
# Utilisons ces indexes triés pour appeler notre GeoDataFrame
# et obtenir les toponymes triés par la taille du cours d'eau :
gdf.loc[lengths.sort_values(ascending=False).index][['TOPONYME']]
Out[19]:
In [20]:
gdf[gdf.TOPONYME == 'Canal du Furon'].index
Out[20]:
In [21]:
geom = gdf.loc[(735, 'geometry')]
In [22]:
print(type(geom))
geom
Out[22]:
In [23]:
print('Longueur: {:.1f}m'.format(geom.length))
Ouverture d'un fichier raster avec rasterio :¶
In [24]:
import rasterio as rio
mnt = rio.open('data/N245E395.tif')
print(mnt.meta)
La projection ne correspond pas, on va reprojeter le fichier raster¶
In [25]:
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.plot import show
dst_crs = 'EPSG:2154' # Lambert 93
transform, width, height = calculate_default_transform(
mnt.crs, dst_crs, mnt.width, mnt.height, *mnt.bounds)
kwargs = mnt.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rio.open('data/N245E395_L93.tif', 'w', **kwargs) as dst:
reproject(
source=rio.band(mnt, 1),
destination=rio.band(dst, 1),
src_transform=mnt.transform,
src_crs=mnt.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest,
)
mnt.close()
In [26]:
mnt = rio.open('data/N245E395_L93.tif')
show(mnt)
Out[26]:
Les deux bibliothèques donnent accès à leurs éléments (géométries d'une part et valeurs du raster d'autre part) sous forme d'objet Python. On va donc facilement pouvoir interagir entre les objets de ces bibliothèques.
Exemple : Création du profil en long d'un cours d'eau à partir du tracé précédement extrait et du mnt reprojeté.
La procédure est relativement simple à expliquer :
- définir une distance entre nos points de relevé (ici on prend arbitrairement 50 points)
- créer des points le long du cours d'eau en utilisant cette distance
- récupérer l'altitude de chacun de ces points
- créer un graphique avec en absisse la distance à la source et en ordonnées l'altitude
In [27]:
# geom est la géométrie récupérée plus tôt
length = geom.length
# combien de points :
nb_points = 50
# quelle distance entre les points :
step = geom.length / nb_points
dists = []
coords = []
for dist in range(0, int(length), int(step)):
pt = geom.interpolate(dist) # type(pt) -> shapely.geometry.Point
# récupération des coordonnées du points sous forme d'un tuple :
c = pt.coords[0]
coords.append(c)
dists.append(dist)
result = []
for i, value in enumerate(mnt.sample(coords)):
result.append((dists[i], value[0]))
On peut maintenant utiliser ces valeurs pour créer un graphique avec matplotlib:
In [28]:
plt.plot([i[0] for i in result], [i[1] for i in result])
plt.title('Profil en long du Canal du Furon')
plt.xlabel('Distance à la source')
plt.ylabel('Altitude')
plt.show()
Out[28]:
Out[28]:
Out[28]:
Out[28]:
Et si on voulait créer une couche de points (où chaque entité serait un de ces points), avec un champ contenant la distance à la source, un champ contenant le toponyme du cours d'eau et un champ contenant l'altitude du point ?
Cette couche va être créée sous forme d'un nouvel objet geopandas et va ensuite être sauvegardée au format GeoJSON.
In [29]:
result = {}
result['geometry'] = []
result['distance_source'] = []
result['altitude'] = []
for dist in range(0, int(length), int(step)):
pt = geom.interpolate(dist)
result['geometry'].append(pt)
result['distance_source'].append(dist)
altitude = list(mnt.sample(pt.coords))[0][0]
result['altitude'].append(altitude)
# Comme avec pandas on peut utiliser un dictionnaire dont les clés sont
# les nom de colonnes et les valeurs associées représentent
# les valeurs des entités de ces colonnes
new_gdf = gpd.GeoDataFrame(result, crs=gdf.crs)
# Toutes nos nouvelles entités vont prendre cette valeur
new_gdf['TOPONYME_COURS_DEAU'] = 'Canal du Furon'
new_gdf.head()
Out[29]:
In [30]:
result = []
for dist in range(0, int(length), int(step)):
pt = geom.interpolate(dist)
result.append({
"geometry": pt,
"distance_source": dist,
"altitude": list(mnt.sample(pt.coords))[0][0],
})
# Cette fois nous avons utilisé une liste
# contenant un dictionnaire pour chaque entité
# ... le résultat est le même qu'avec la méthode prédécente
new_gdf = gpd.GeoDataFrame(result, crs=gdf.crs)
# Toutes nos nouvelles entités vont prendre cette valeur
new_gdf['TOPONYME_COURS_DEAU'] = 'Canal du Furon'
new_gdf.head()
Out[30]:
In [31]:
new_gdf.plot()
Out[31]:
Enregistrement du résultat¶
Transformation des géometries en coordonnées géographiques (EPSG:4326) avant d'enregistrer le fichier au format GeoJSON :
In [33]:
new_gdf.to_crs(epsg=4326).to_file('data/points_canal_furon.geojson', driver='GeoJSON')
In [34]:
%cat data/points_canal_furon.geojson
Plus de manipulations géospatiales avec geopandas¶
On a vu que geopandas permettait de manipuler les données sous la même forme que pandas (lignes / colonnes). Cette possibilité va pouvoir être utilisée pour appliquer des opérations à l'ensemble de la colonne geometry.
On peut par exemple créer un buffer autour de chaque point de la manière suivante :
In [35]:
# On peut appeler la méthode sur l'ensemble de l'objet ou seulement sur sa colonne geometry, le résultat est le même
# Dans les deux cas c'est une GeoSeries (une colonne et un index) qui est renvoyée :
new_gdf.buffer(100).head()
print('new_gdf.buffer(100) == new_gdf.geometry.buffer(100) -> {}'
.format(np.all(new_gdf.buffer(10) == new_gdf.geometry.buffer(10))))
Out[35]:
On peut utiliser cette colonne pour remplacer la colonne de géométries qui contient nos points :
In [36]:
new_gdf.geometry = new_gdf.geometry.buffer(100)
In [37]:
new_gdf.head()
new_gdf.tail()
Out[37]:
Out[37]:
Plus de manipulations géospatiales avec geopandas (suite)¶
In [38]:
# Ouverture de la couche des communes de la région Auvergne-Rhône-Alpes
comm = gpd.read_file('data/COMMUNE_auvergne_rhone_alpes.shp')
# Le système de coordonnées est le même que nos cours d'eau
assert comm.crs == gdf.crs
In [39]:
# Aggrégation / fusion des communes en départements
# On choisit les colonnes à conserver, la colonne à utiliser
# pour l'aggrégation, et une fonction d'aggrégation si besoin
dep = comm[['geometry', 'NOM_DEPT', 'SUPERFICIE', 'POPULATION']].dissolve(by='NOM_DEPT', aggfunc='sum')
dep.reset_index(drop=False, inplace=True)
# Affichage du résultat
dep.head(4)
dep.plot()
Out[39]:
Out[39]:
In [40]:
geom_furon = gdf.loc[gdf.TOPONYME == 'Canal du Furon', 'geometry'].values[0]
# Quelles sont les communes traversées par le Canal du Furon ?
comm[comm.intersects(geom_furon.buffer(1000))]
Out[40]:
In [41]:
# Calculer la densité de population sur la couche des communes :
comm['DENSITE_POPULATION'] = comm['POPULATION'] / comm['SUPERFICIE']
comm.plot('DENSITE_POPULATION', scheme='Quantiles', k=8, figsize=(16,16), legend=True)
Out[41]:
Attention geopandas ne dispose toutefois pas de 'vrais' fonctionnalités 'carto' et il existe des bibliothèques qui offrent une API dédiée à la cartographie en Python.¶
Sans se plonger dans la cartographie thématique, on note qu'il est facilement possible des faire des cartes choroplètes (avec geopandas, ou folium ensuite par exemple) ou des cartes de symboles proportionnels avec quelques opérations supplémentaires.¶
Affichage de nos 50 buffer sur un fond Leaflet :¶
In [42]:
import folium
from pyproj import Proj, transform
# On va récupérer les coordonnées du centroide du cours d'eau pour centrer la carte
# (on doit transformer ces coordonnées du Lambert 93 vers des coordonnées en latitude-longitude)
c = geom.centroid.xy
x,y = transform(Proj(init="epsg:2154"), Proj(init="epsg:4326"), c[0][0], c[1][0])
# On instancie un objet Map avec le fond de notre choix
# et on lui ajoute direction nos entitées polygonales:
map_example = folium.Map(
location=(y, x),
tiles = "Stamen Toner",
zoom_start = 11).add_child(
folium.features.GeoJson(new_gdf.to_crs(epsg='4326').to_json()))
In [43]:
# Export du document HTML
map_example.save('tmp/plot_data.html')
# Affichage dans IPython
IFrame('tmp/plot_data.html', width="700px", height="400px")
Out[43]:
Dans la partie suivante nous verrons comment créer une page réactive utilisant des données coté serveur ainsi qu'une solution rapide pour exposer une API réalisant des calculs personnalisés par exemple.
Bien que ça ne soit pas le cas dans les exemples suivants, il est possible d'intégrer des données géo-spatiales dans développements web comme ceux qui vont suivre.