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 de QGIS

Données raster

Carte de probabilité de pluie

Données 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)

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 de l'information à partir de données vectorielles

rasterstats : statistiques sommaires de jeux de données raster géospatiaux selon des géométries vectorielles (https://github.com/perrygeo/python-rasterstats)

Extraire des valeurs raster avec rasterstats

Pour les vecteurs de points :

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

Pour les vecteurs de 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