Travailler avec des données raster

Travailler avec des données géospatiales en Python

Joris Van den Bossche

Open source software developer and teacher, GeoPandas maintainer

Illustration d'une grille de valeurs

Source de l'image : documentation QGIS

Données raster

Carte de probabilité de pluie

Raster à bandes multiples

Le paquet rasterio

import rasterio

Liaisons « pythoniques » à GDAL
Lecture et écriture de fichiers raster
Outils de traitement (masquage, reprojection, rééchantillonnage, …)

https://rasterio.readthedocs.io/en/latest/

Ouvrir un fichier raster

import rasterio

src = rasterio.open("DEM_world.tif")

Métadonnées :

src.count

src.width, src.height

(4320, 2160)

Les données raster = tableau numpy

array = src.read()

Tableau numpy standard :

array

array([[[-4290, -4290, -4290, ..., -4290, -4290, -4290],
        [-4278, -4278, -4278, ..., -4278, -4278, -4278],
        [-4269, -4269, -4269, ..., -4269, -4269, -4269],
        ...,
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804]]], dtype=int16)

Tracer un jeu de données raster

Avec la méthode rasterio.plot.show() :

import rasterio.plot

rasterio.plot.show(src, cmap='terrain')

Extraire des informations depuis des données vectorielles

rasterstats : statistiques synthétiques de rasters géospatiaux basées sur des géométries vectorielles (https://github.com/perrygeo/python-rasterstats)

Extraire des valeurs raster avec rasterstats

Pour des points :

rasterstats.point_query(geometries, "path/to/raster", 
                        interpolation='nearest'|'bilinear')

Pour des polygones :

rasterstats.zonal_stats(geometries, "path/to/raster",
                        stats=['min', 'mean', 'max'])

Extraire des valeurs raster avec rasterstats

result = rasterstats.zonal_stats(countries.geometry, "DEM_gworld.tif", 
                                 stats=['mean'])

countries['mean_elevation'] = pd.DataFrame(result)

countries.sort_values('mean_elevation', ascending=False).head()

            name    continent                     geometry  mean_elevation
157   Tajikistan         Asia  POLYGON ((74.98 37.41, ...      3103.231105
85    Kyrgyzstan         Asia  POLYGON ((80.25 42.34, ...      2867.717142
24        Bhutan         Asia  POLYGON ((91.69 27.77, ...      2573.559846
119        Nepal         Asia  POLYGON ((81.11 30.18, ...      2408.907816
6     Antarctica   Antarctica  (POLYGON ((-59.57 -80.04...     2374.075028
..           ...          ...                          ...             ...

Passons à la pratique !

Travailler avec des données géospatiales en Python