Wiki du Master Géographies Numériques

L'émergence des SIG a permis aux ingénieurs, aux géographes, aux biologistes, et autres professionnels, de travailler en couches, c'est-à-dire de travailler avec la superposition de plusieurs informations relatives à un lieu (densité démographique, utilisation du sol, topologie, et toute autre caractéristique ayant un x et un y) et de réaliser des analyses à partir de la relation entre elles. La programmation permet de réduire le temps nécessaire à la réalisation de ces tâches et d'augmenter l'échelle spatiale des lieux où ces analyses sont effectuées à partir d'automatisations que l'homme ne peut pas réaliser seul.

L'automatisation proposée par la programmation peut être réalisée, notamment dans les logiciels SIG traditionnels, aussi bien sous forme de modèles graphiques que sous forme écrite. Les modèles graphiques représentent une manière plus visuelle d'appliquer l'automatisation, parmi lesquels nous pouvons citer ModelBuilder d'Esri et le Modeleur Graphique de Qgis. A l'écrit, nous pouvons voir l'adaptation de langages traditionnels pour le monde des SIG comme Arcpy et PyQgis, en plus ces deux logiciels insèrent dans leurs interfaces l'option d'ouvrir une console propre pour appliquer ces langages.

En plus de l'insertion de Python dans les logiciels, il est également possible de constater la création de librairies permettant l'utilisation de couches vectorielles et même de rasters. Ces bibliothèques permettent non seulement la visualisation de ces fichiers, mais aussi l'application d'outils de géotraitement et la production de cartes.

Parmi les exemples de librairies, citons Geopandas, qui est une extension du Pandas avec l'ajout de la composante géométrique des données et peut donc ouvrir des fichiers shapefile et geopackage, par exemple. Une autre librairie bien connue est GDAL/OGR, qui permet de lire non seulement les fichiers vectoriels mais aussi les rasters. Enfin, nous pouvons également mentionner la librairie rasterio, axée sur la manipulation des fichiers rasters. De tels outils permettent la réalisation de fonctions connues comme Clip, Buffer, Merge, Centroid, Dissolve, etc. En outre, ils permettent également d'effectuer des calculs tels que la surface, la longueur, la distance et autres. Toutes ces librairies disposent d'une documentation expliquant comment les utiliser.

On peut donc observer que la programmation fait de plus en plus partie du monde de SIG, soit comme outil auxiliaire, soit comme outil SIG lui-même. Voici quelques exemples de comment utiliser la programmation et l’automatisation dans le monde de SIG.

Pour accéder au ModelBuilder il faut cliquer sur l’onglet « Analysis » et chercher dans l’espace « Geoprocessing » son icon. Une fois cliqué sur l’icon, un nouvel onglet apparaitra.

Source : ESRI.

Dans le nouvel onglet, on trouve plusieurs commandes qu’on peut utiliser avec ModelBuilder, cependant la plus important est l’espace « Insert » qui contient les Itérateurs, les Utilities et les opérateurs de Logique.

Source : ESRI.

Les noms des Itérateurs dans ce tuto seront en anglais, justement pour faciliter la recherche sur internet quand il y a des doutes. Parce qu’il y a plus de sources d’aide en anglais qu’en français. Cependant, la liste de itérateurs avec sa traduction est montrée ci-suivante.

1. For est Pour ;
2. While est Pendant ;
3. Iterate Feature Selection est Itérer la sélection d’entités ;
4. Iterate Row Selection est Itérer la sélection de lignes ;
5. Iterate Fields est Itérer les champs ;
6. Iterate Field Values est Itérer les valeurs de champs ;
7. Iterate Multivalue est Itérer les valeurs multiples ;
8. Iterate Dataset est Itérer les jeux de données ;
9. Iterate Features Classes est Itérer les classes d’entités ;
10. Iterate Files est Itérer les fichiers ;
11. Iterate Layers Itérer sur les couches ;
12. Iterate Rasters est Itérer les rasters ;
13. Iterate Tables est Itérer les tables ;
14. Iterate Workspaces est Itérer les espaces de travail.

Comme Itérateurs on a beaucoup d’options et il est possible d’utiliser des itérations plus connues comme For et While mais aussi des itérations liées aux entités du monde SIG. Il est possible de réaliser une itération par rapport à table attributaire, on utilise « Iterate Feature Selection » pour les couches ou entités et « Iterate Row Selection » pour les tables. For est intéressant pour des analyses répétitives avec une limite, comme par exemple dans l’image suivante, on réalise plusieurs zones tampon de 500 à 2000.

Source : (ESRI).

On peut aussi réaliser des itérations par rapport aux champs d’une couche présents dans l’espace Contenu avec l’itérateur « Iterate Fields », ça permet de les filtrer par type (Texte, Réel, Court, Long, étc) et aussi de les sélectionner par nom. À travers l’itérateur « Iterate Fields Values », on peut aller plus loin et réaliser des itérations à propos des valeurs de chaque champ. L’itérateur « Iterate Multivalue » est capable d’itérer des plusieurs classes d’entités dedans différents fichiers geodatabase. Cet itérateur va parcourir chaque classe d’entité et réaliser une action sur elle. Dans l’image suivante on peut voir que chaque classe d’entité (InputFC1, InputFC2, InputFC3) sera projetée à travers l’outil « Project ».

Source : (ESRI).

L’itérateur « Iterate Dataset » permet de travailler avec des dossiers ou fichiers geodatabase. Ainsi, c’est possible d’itérer sur les fichiers contenus et aussi de les filtrer par type de fichier (Feature, Raster, TIN, étc).

Ainsi comme l’itérateur précèdent, l’itérateur « Iterate Feature Classes » il est possible de réaliser des itérations au niveau des classes d’entités dans un jeu de classes d’entités (soit un fichier geodatabase soit un dossier normal contenant des shapefiles). Cet itérateur permet aussi de réaliser un filtre par type des d’entités, c’est-à-dire qu’on peut choisir séparément des points, lignes ou polygones, à travers l’espace « Feature Type ». En plus, on peut aussi réaliser un filtre par le nom des fichiers à travers l’espace « Wildcard » avec l’aide du caractère spécial « * », regardez quelques exemples :

1. « Grand* » - on filtre tous les entités qui commencent par « Grand » ;
2. « *Grand* » - on filtre tous les entités qui ont le mot « Grand » au milieu de son nom ;
3. « *Grand » - on filtre tous les entités qui finissent par « Grand », sans considérer l’extension du fichier.

Dans l’image suivante on peut trouver un exemple où l’itérateur « Iterate Feature Classes » sélectionne les entités du type polygone d’un geodatabase appelé « Output.gdb ». Ensuite, pour chaque entité il est ajouté et calculé un nouveau champ à travers l’outil « Add Field » et « Calculate field », respectivement.

Source : (ESRI).

L’itérateur « Iterate Files » permet de travailler avec des fichiers dans un dossier. Également comme « Iterate Feature Classes », il est possible de réaliser des filtres à travers des mots (utilisation d’un mot plus le caractère « * ») et aussi l’extension du fichier (par exemple, pour travailler avec des fichiers du type texte, il faut remplir « TXT »).

Source : (ESRI).

On peut aussi réaliser des itérations des couches dans la Carte d’un fichier *.aprx, pour le faire on utilise l’itérateur « Iterate Layers ». Cet itérateur permet de réaliser des filtres au niveau des noms des couches (Wildcard), par rapport au type de couche (Layer type) et aussi par rapport si les couches sont couchées ou pas. En choisissant l’option Feature Layer comme type, l’itérateur augmente les options des filtres. On peut filtrer la source de la couche (Workspace Type) et le type (Feature Type, par exemple, points, lignes ou polygones).

Les trois derniers itérateurs fonctionnent de façon similaire, il faut leur donner un dossier ou un fichier geodatabase sur l’espace Workspace et l’itérateur va parcourir chaque type de fichier. L’itérateur « Iterate Rasters » va parcourir les fichiers du type raster et on peut filtrer par nom (Wildcard) et aussi par type de raster (Raster Format). L’itérateur « Iterate Tables » va parcourir les tables et il permet aussi de filtrer para nom et type de table. L’itérateur « Iterate Workspaces » permet de parcourir des dossiers dedans un dossier fourni et comme les autres il est possible de filtrer par nom et type de dossier.

La Bibliographie utilisée par ce tutoriel est la page d'ESRI Pour la consulter il faut utiliser les liens suivants :

1. Pour consulter les itérateurs : https://pro.arcgis.com/fr/pro-app/latest/tool-reference/modelbuilder-toolbox/examples-of-using-iterators-in-modelbuilder.htm#GUID-C0E9F998-C1DE-4D52-9FCE-13543FBC0FCB ;
2. Pour consulter tous les outils possibles sur ModelBuilder : https://pro.arcgis.com/fr/pro-app/latest/tool-reference/modelbuilder-toolbox/an-overview-of-the-modelbuilder-toolbox.htm ;

La première chose à faire avant d’utiliser les commandes de la librairie GeoPandas est l’installation et pour la faire il faut savoir quelle plateforme IDE vous utilisez (Jupyter Notebook, Google Colab, Anaconda, etc). Dans tout le cas la page de documentation de GeoPandas fourni des consignes comment faire son installation : https://geopandas.org/en/stable/getting_started/install.html. Pour Jupyter Notebook il faut faire comme dans l’image suivante et installer tous les libraires nécessaires avant d’installer geopandas, c’est-à-dire shapely, fiona, pyproj et packaging.

Ensuite on peut installer geopandas. Les librairies mapclassify et folium seront utile pour la visualisation des couches avec un fond de carte.

Une fois installé, il faut importer cette librairie en écrivant « import geopandas ». Ensuite on peut finalement commencer à utiliser les commandes du GeoPandas. Cette librairie permet de travailler avec les fichiers du type GeoPackage, GeoJSON et Shapefile et tous les autres types contenant une table attributaire et une geometrie. Dans l’image on voit que d’autres librairies ont étés installés, c’est normal parce qu’on aura besoin pour la séquence du tutoriel.

Pour lire des fichiers il faut utiliser la commande « geopandas.read_file() » et y ajouter entre les parenthèses le chemin du fichier avec son nom et son extension (*.shp, *.geojson ou *.gpkg par exemple). Au final on écrit geopandas.read_file(« chemin du fichier/nom.extension »). On peut aussi stocker dans une variable pour faciliter l’utilisation de cette commande, « variable = geopandas.read_file(« chemin du fichier/nom.extension »). Le résultat est une table avec les champs de la table attributaire comme colonne et à la fin une dernière colonne qui montre la géométrie (polygon, ligne ou point par exemple) comme dans l’image ci-suivante.

GeoPandas permet aussi de visualiser les données, on doit utiliser la commande plot(). En donnant le nom d’une colonne, on peut avoir une symbologie par valeur unique et on a aussi la possibilité de montrer la légende, comme dans l’image ci-dessous :

On peut aussi utiliser un fond de carte, mais on doit utiliser la commande .explore(). Cependant, il faut installer et importer les librairies folium, matplotlib et mapclassify, comme on avait dit avant. Le résultat est une carte interactive avec la couche sur le fond de carte du OpenStreetMap. Avec cette carte, c’est possible de zoomer et bouger l’écran.

La commande to_file() sert à deux possibilités, soit créer des nouveaux fichiers soit convertir des fichiers existants. Il faut seulement donner comme arguments le chemin, avec le nom et extension, et l’extension (argument « driver »). Un exemple est montré dans l’image suivante, dans ce cas on faire une conversion Shapefile > Geopackage. Pour quelques formats il faut savoir la forme abregée qui GeoPandas accepte, ce que peut être consulter dans la page documentation (Lien : https://geopandas.org/en/stable/docs/user_guide/io.html ).

En plus de lire, écrire et convertir des fichiers, la librairie GeoPandas contient plusieurs fonctions géospatiales. Voici une liste de fonctions géospatiales : .area, .boundary(), .centroid(), .distance(), to_crs(), contains(), intersects(), disjoint(), etc.

La fonction to_crs() permet de convertir le système des coordonnées du fichier. S’on ne connait pas le système de coordonnées de la couche, on peut utiliser la commande .crs et on aura tous les informations par rapport à ces paramètres, comme on peut voir dans l’image suivante.

Ensuite, on peut convertir de WGS84 à Lambert 93 avec .to_crs() en donnant le code EPSG correct. On peut stocker la couche nouvelle dans une autre variable.

Comme GeoPandas est une sous-classe de Pandas on peut utiliser toutes les fonctions Pandas. Ça nous donne la possibilité de faire des sélections par attribut en utilisant les commandes .loc[] e .iloc[]. L’image suivante contient un exemple de comment sélectionner les entités avec une surface plus grand que 5000 km².

On peut faire aussi des filtres avec des champs du type string et sélectionner seulement les départements avec «-» dans son nom. La fonction « .str » sert à convertir les valeurs de la colonne « NOM_M » en type STRING ou Texte en français.

D’autres fonctions et commandes sont possible avec la librairie GeoPandas et on peut les trouver sur son page de Documentation (Lien : https://geopandas.org/en/stable/docs.html). Dans cette page, il contient les informations utilisées pour réaliser ce tutoriel.

A travers la librairie GDAL/OGR on peut travailler avec des fichiers vectoriels. Il faut utiliser OGR. Avant de commencer il faut dire que la librairie OGR n’est pas si « user friendly » que GeoPandas et elle travaille d'une manière qui demande à l'utilisateur de fournir toutes les informations, alors que GeoPandas est capable de comprendre tout seul.

L’importation est assez simple.

Après l’importation on peut lire des fichiers vectoriel à partir de la fonction « ogr.Open(chemin du fichier). On peut stocker l’ouverture dans une variable. Ainsi comme GeoPandas, OGR est capable de travailler avec beaucoup type de fichier vectoriel, la liste de tous les fichiers est trouvée dans le lien suivant https://gdal.org/drivers/vector/index.html.

Pour commencer à travailler plus à fond avec cette librairie, et c’est pour cette raison qu’elle est un peu plus difficile que GeoPandas, il faut définir ce qui est une « layer ». Pour faire ça il faut utiliser la fonction « .GetLayer() ».

A partir de « layer » on peut avoir des informations par rapport au fichier vectoriel, comme le système des coordonnées. En utilisant la fonction « .GetSpatialRef() » on obtient le système des coordonnées et en appliquant « .ExportToWkt() » on obtient l’info du système des coordonnées dans une façon plus facile de lire.

On peut aussi extraire la géométrie du fichier vectoriel à travers de « layer ». Cependant il faut une étape avant, c’est la création de la « feature » de « layer », on utilise la fonction « GetFeature() ». Ensuite on utilise la fonction « GetGeometryRef() » et la fonction « ExportToWkt() » pour avoir la géométrie. On stocke la géométrie dans la variable « geom ».

Avec la géométrie, on peut appliquer des calculs du périmètre ou de la surface du fichier vectoriel. Le périmètre est obtenu à travers deux fonctions : « .Boundary() » et « Length() ». Comme le système de coordonnées est Lambert 93, l’unité est mètre.

La surface est déterminée à partir de la fonction « Area() ». Le résultat est donné en fonction du système de coordonnées, comme il s’agit de Lambert 93, l’unité est m².

On peut aussi réaliser des buffers avec la géométrie en utilisant l’outil « .Buffer() ».

Il y a plusieurs outils dans la librairie OGR, cependant elle demande un niveau très haut de connaissance. Pour avoir acquérir plus de connaissance par rapport au module osgeo.ogr spécialement, il est possible de consulter sa page de documentation https://gdal.org/api/python/osgeo.ogr.html#module-osgeo.ogr.

La librairie GDAL/OGR permet de travailler non seulement avec des fichiers vectoriels mais aussi des fichiers rasters. Pour les fichiers rasters on utilise GDAL. Pour l’utiliser il faut l’importer et aussi la librairie matplotlib pour visualiser les fichiers.

Pour ouvrir un fichier raster il faut utiliser la fonction « gdal.open() » et la donner le chemin du fichier raster. GDAL est capable de lire plusieurs types de fichiers raster, la liste contenant les types est trouvé dans le lien suivant https://gdal.org/drivers/raster/index.html. On va utiliser un fichier du type BD Alti IGN réfèrent au autour de Saint-Étienne.

On peut obtenir des information du raster à partir de la fonction « .GetGeoTransform() ». Elle nous fournit comme résultat contenant comme premier et quatrième élément, respectivement, les coordonnées x et y du système des coordonnées référents au point d’origine à gauche et en haut, elle fournit aussi la résolution horizontal (deuxième élément) et la résolution vertical (dernier élément, et en négatif par defaut).

Le système des coordonnées du fichier raster peut être obtenu avec l’utilisation de la fonction « .GetProjection() ». Le résultat contient plusieurs informations comme le nom, le code EPSG, l’unité, etc.

Avec GDAL, on peut réaliser aussi des analyses en utilisant un MNT comme Slope, Hillshade et Aspect. Pour le faire il faut travailler avec l’outil « osgeo.gdal.DEMProcessing(destName, srcDS, processing, computeEdge) ». L’argument « destName » est le chemin du fichier résultante, l’argument « srcDS » est le fichier original et l’argument « processing » est le type d’analyse que peut être : “hillshade”, “slope”, “aspect”, “color-relief”, “TRI”, “TPI” et “Roughness”. L’argument « computeEdge » est pour définir si l’analyse sera réalisée aux niveaux des limites du fichier raster. Pour avoir plus d’explication on peut consulter la page de documentation de l’outil sur https://gdal.org/api/python/osgeo.gdal.html.

On va faire trois exemples, le slope, l’hillshade et l’aspect. Ensuite on va les visualiser en utilisant la commande plot de matplotlib (« figure » pour créer les axis, « imshow » pour définir le fichier de l’analyse réalisée à montrer, « colorbar » pour montre la légende par rapport à la couleur du raster et « show » pour montrer le résultat).

Pour le slope il faut écrire « osgeo.gdal.DEMProcessing('slope_sainte.tif', sainte_bd_alti, “slope”, computeEdges = True) » et on stocke sur la variable « slope ». À travers cette ligne de code on va créer un fichier « slope_sainte.tif » dans le même chemin du fichier Jupyter Notebook contenant les valeurs des pentes.

Pour le hillshade, le méthode est similaire et il faut écrire « osgeo.gdal.DEMProcessing('hillshade_sainte.tif', sainte_bd_alti, “hillshade”, computeEdges = True) » et stocker dans la variable « hillshade ». On crée un fichier qui s’appelle « hillshade_sainte.tif ».

Pour l’aspect, on fait pareillement et on écrit « osgeo.gdal.DEMProcessing('aspect_sainte.tif', sainte_bd_alti, “aspect”, computeEdges = True) » et on stocke sur la variable « aspect ». Comme dans les dernières analyses, on cree aussi un fichier qui s’appelle « aspect_sainte.tif ».

Une fois crées les fichiers des analyses comme slope, hillshade ou aspect, on peut les ouvrir sur QGIS ou ArcGIS PRO. Cependant il faut fermer tous les fenêtres Jupyter Notebook, puisque les fichiers ne sont plus connectés à la plateforme python.

La librairie GDAL/OGR contient plusieurs fonctions qui peuvent être utilisé avec fichier vectoriels et rasters, et pour les découvrir il faut consulter la page de documentation sur le lien suivant https://gdal.org/index.html. Et pour les fichiers raster, on peut consulter la page https://gdal.org/programs/index.html. Dans cette page, il contient les informations utilisées pour réaliser ce tutoriel.

Comme pour la librairie GeoPandas, il faut installer la librairie Rasterio. Pour en faire on doit écrire seulement « !pip install rasterio -q ». Pour quelques problèmes d’installation qui peuvent apparaître, Rasterio fournit une page en montrant comme installer proprement https://rasterio.readthedocs.io/en/stable/installation.html#.

Ensuite, on installe la librairie à travers la commande import, comme dans l’image suivante. On importe aussi la fonction plot de Rasterio pour nous aider dans la visualisation des rasters.

Après l’installation et l’importation, on peut commencer à utiliser Rasterio. Pour ouvrir des fichiers rasters il faut utiliser la fonction « rasterio.open » et donner le chemin du fichier. Et comme dans GeoPandas, on peut stocker l’ouverture d’un fichier dans une variable (« image », dans ce cas). La fonction rasterio.open permet de lire rasters sous le format Tiff, GeoTiff, ASCII grid et Erdas Imagine par exemple. Comme cette librairie est basée sur GDAL, elle peut lire tous les fichiers possible sur GDAL (https://gdal.org/drivers/raster/index.html).

Pour savoir les propriétés du raster, la fonction « .meta » est la solution. A travers d’elle on a accès aux informations metadata du fichier comme : driver (format du raster), dtype (type de espace du raster), nodata (présence de nodata), width (largeur du raster), height (hauteur du raster), count (nombre de bands), crs (info par rapport au système des coordonnées), transform (matrice de transformation affine qui correspondre les emplacement des pixels en coordonnées, en col et ligne, à des positions spatiales, x et y). L’image utilisé est une composition vrai couleur de l’autour de Lyon, c’est pour cette raison qui dans « count » il y a le valeur 3 (3 bands : rouge, vert et bleu).

La visualisation de l’image fonctionne en utilisant la fonction « plot.show(variable qui stock l’ouverture) ». On peut visualiser l’image dans une graphe qui utilise les coordonnées longitude et latitude comme axis x et y, respectivement.

Comme exemple d’application on va voir comment réaliser une extraction par masque. Pour ce but il faut réaliser trois étapes avant. La première est l’importation du GeoPandas pour nous aider avec l’ouverture des fichiers vectoriel. Ensuite, on importe aussi la librairie json pour extraire les coordonnées du fichier vectoriel. Et on importe la fonction rasterio.mask.mask de rasterio.

On utilise geopandas.read_file() pour ajouter la couche de la ville de Lyon. Cette couche est en format *.shp. On va convertir le format en exportant sous le format *.geojson et ensuite on va stocker la couche résultante dans la variable « couper » parce c’est elle qui servira comme masque. Ensuite il faut qu’on extraie les coordonnées du fichier geojson, on utilise la fonction json.loads() et on défine que les coordonnées sont dedans « features » et « geometry » dans la structure du fichier geojson. Les coordonnées seront stockées dans la variable « coords ».

Comme on peut voir dans l’image suivante les coordonnées ont été bien stockées dans la variable « coords ». Cependant, il faut remarquer qu’elles sont encore dans une structure du format json.

Avec toutes les préparations finies, on peut finalement appliquer la fonction « mask » de rasterio. Il faut qu’on donne comme arguments la variable qui stocke le raster (image), la variable contenant les coordonnées (coords) et qu’on écrire « crops = True » (cet arguments permet de couper l’image). On va stocker le resultat dans deux variables : masque_img qui va stocker l’image coupé et masque_transform qui va strocker les information par rapport à transform (une information de metadata).

Ensuite, on va mettre à jour les informations de metadata qu’on a copié dans la variable « masque_meta ». Pour finir, on enregistre le fichier resultant, comme on a donné seulement le nom du fichier, il sera enregistré dans le même chemin du fichier Jupyter Notebook. Le « w » sert à définir qu’il s’agit d’une création de fichier et il vient de la mot anglais « write ».

On ouvre le fichier avec la fonction « rasterio.open() » et ensuite on utilise « plot » pour visualiser le résultat final.

On vérifie aussi les informations metadata à travers la commande « .meta ». On peut voir que les informations ont été bien mis à jour.

La librairie Rasterio permet de réaliser diverses fonctions et ainsi comme la librairie GeoPandas, elle a une page de documentation contenant tous ses fonctions et aussi l’explication des arguments nécessaires. Le lien de la page est le suivant https://rasterio.readthedocs.io/en/stable/api/rasterio.html. Dans cette page, il contient les informations utilisées pour réaliser ce tutoriel.

Matheus Gimenez Saraiva - 2022/2023