Analyses géospatiales avec le logiciel commercial ArcGIS d’ESRI et le paquet arcpy. Ce portable ne peut pas être exécuté dans Jupyter et les blocs de code nécessitent la bibliothèque commerciale arcpy d’ESRI.
Présentation¶
Le logiciel d’ESRI ArcGIS Pro est livré avec son propre environnement conda qui permet de travailler avec la bibliothèque commerciale arcpy. arcpy est accessible soit directement dans ArcGIS Pro ou en utilisant l’environnement conda de l’ESRI, qui peut être géré par ArcGIS Pro. Puisque arcpy est lié à l’utilisation d’un environnement commercial, il ne peut être consulté que dans des cahiers de jupyter directement hébergés à esri.com (mais ceux-ci ne fonctionneront pas avec les exemples de cette section). Cette page explique l’utilisation de base de arcpy à l’externe IDEs (p. ex., PyCharm, VS Code, ou Spyder), l’intégration des licences, et la fonctionnalité de base de arcpy.
Historique¶
Pourquoi travailler avec arcpy dans Python reste-t-il un outil puissant même s’il existe de nombreuses restrictions de licence et de plateforme ?
Si vous travaillez pour une entreprise avec des services d’ingénierie, l’entreprise utilise probablement ArcGIS en raison de sa popularité et service de soutien commercial. Dans ce cas, une entreprise (ou un établissement de recherche) couvre les frais de licence et il est tout aussi probable qu’un système d’exploitation Windows soit utilisé pour des raisons similaires. Parmi les fonctionnalités populaires du logiciel ArcGIS figurent Raster Calculator, Shapefile Feature, ou des outils pour l’analyse statistique des bases de données géographiques. Grâce à arcpy, ces outils ArcGIS populaires peuvent être intégrés dans les scripts Python, ce qui permet d’automatiser les flux de travail et d’accroître considérablement l’efficacité.
Cette page montre les bases pour manipuler les données de raster et de shapefile dans Python avec arcpy. Il y a beaucoup plus de méthodes implémentées à arcpy et seuls les fondamentaux sont expliqués ici.
Utiliser arcpy avec les IDE externes¶
Interprète de configuration¶
La section d’installation Python explique how to set up PyCharm IDE avec un environnement conda. Pour créer un nouvel environnement de conda dans PyCharm avec l’environnement de conda d’ESRI, la Location doit être définie différemment (voir la screenshot]) :
Si ArcGIS Pro a été installé dans le monde entier par l’administrateur du système, utiliser:
%PROGRAMFILES%\ArcGIS\Pro\bin\Python\Scripts\propySi ArcGIS Pro a été installé pour les utilisateurs individuels, utiliser:
%LOCALAPPDATA%\Programs\ArcGIS\Pro\bin\Python\Scripts\propy
Il se peut qu’il ne soit pas nécessaire d’ajouter le répertoire propy. En outre, trouvez le conda.exe ou python.exe dans les répertoires ci-dessus et définissez-le dans le champ Conda exécutable. Finaliser la création du nouvel environnement en cliquant sur Créer.
Pour installer d’autres paquets, suivez les descriptions fournies par Esri.
Importer arcpy et ses modules¶
arcpy est livré avec ses propres classes pour les objets de données géospatiales (par exemple, arcpy.Raster pour les données maillées) et les modules pour la cartographie (arcpy.mp), l’émulation de l’analyste Spatial (arcpy.sa), ou l’accès aux données (arcpy.da). Une liste complète est fournie sur le site Web du développeur. Cette page comporte des blocs de code à l’aide de arcpy et arcpy.sa, où les objets Spatial Analyst sont importés à l’aide de * pour permettre un accès direct, par exemple, à Con() (au lieu de arcpy.sa.Con()), ce qui correspond à la déclaration conditionnelle de arcpy pour les rasters.
import arcpy
from arcpy.sa import *Configuration arcpy Espace de travail¶
Les calculs géospatials peuvent produire de nombreux produits secondaires, qui peuvent être lourds de taille. Pour mieux contrôler quelles données sont générées par arcpy et où il est bon de définir un espace de travail pour chaque script arcpy avec :
arcpy.env.workspace = "C:\\workspace\\" # or use os.path.dirname(__file__) to go to the script directoryIl peut également être utile de permettre l’écrasement de fichiers déjà existants avec le même nom en raison de la taille du fichier. Ce comportement peut être souhaité ou non et peut être contrôlé comme suit.
arcpy.gp.overwriteOutput = True # enable overwriting
arcpy.gp.overwriteOutput = False # disable overwritingDéfinir les étendues spatiales¶
Les ensembles de données géospatiales peuvent avoir une grande étendue sans qu’aucune donnée ne soit écrite en grande partie. Le traitement des cellules de données vides peut entraîner un temps de calcul inutilement long. Par conséquent, il est conseillé de limiter l’étendue du calcul avec:
arcpy.env.extent = "MAXOF" # uses the combined extent of all input datasets
arcpy.env.extent = "MINOF" # uses only the overlap of all input datasets
arcpy.env.extent = arcpy.Extent(arcpy.Raster("base.tif")) # uses the controlled extent of a raster
arcpy.env.extent = "Xmin YMin XMax Ymax" # imposes user-defined minimum and maximum coordinatesLicences de vérification¶
Plusieurs méthodes arcpy ont besoin de licences comme Spatial Analyst ou 3D. Dans les scripts autonomes, les licences peuvent être activées (checked out) avec arcpy.CheckOutExtension('NAME') et désactivées (checked in) avec arcpy.CheckInExtension('NAME'). Compte tenu de l’orientation de l’objet, les opérations arcpy doivent être intégrées dans les fonctions ou les méthodes des classes. Par conséquent, il est recommandé d’envelopper les fonctions ou les méthodes en utilisant arcpy avec decorators qui activent les licences nécessaires. Le bloc de code suivant fournit un wrapper pour activer une licence Spatial Analyst pour une fonction.
def spatial_license(func):
def wrapper(*args, **kwargs):
arcpy.CheckOutExtension('Spatial')
result = func(*args, **kwargs)
arcpy.CheckInExtension('Spatial')
return result
return wrapperErreurs de suivi¶
La section Python Errors, Logging, and Debugging fournit des instructions utiles pour le dépannage des erreurs dans l’utilisation du code. Pour identifier les problèmes liés aux scripts arcpy, il est recommandé d’utiliser une fonction d’emballage supplémentaire, qui écrit des erreurs arcpy dans un fichier journal (logger). Les expressions logger.*(MESSAGE) peuvent également être remplacées par print(MESSAGE).
import logging
def err_info(func):
def wrapper(*args, **kwargs):
arcpy.gp.overwriteOutput = True
logger = logging.getLogger("logfile")
try:
return func(*args, **kwargs)
except arcpy.ExecuteError:
logger.info(arcpy.GetMessages(2))
arcpy.AddError(arcpy.GetMessages(2))
except Exception as e:
logger.info(e.args[0])
arcpy.AddError(e.args[0])
except:
logger.info(arcpy.GetMessages())
return wrapperAnalyste spatial et opérations de restauration¶
Principes de base¶
arcpy offre différentes options pour charger géospatial rasters (données maillées), qui peuvent avoir différents formats tels que Esri Grid (sans fin, le raster est un dossier avec d’autres fichiers), GeoTIFF, DAT et bien d’autres. Le script suivant charge une profondeur de flux Grille raster h et une vitesse de flux Grille raster u:
h = arcpy.Raster("geodata/input/rasters/h")
u = arcpy.Raster("geodata/input/rasters/u")Le Froude number peut être calculé pixel par pixel à partir de la profondeur et de la vitesse du flux et de la constante de gravité =9.81 m/s. Le script suivant calcule le nombre de Froude pour tous les pixels où la profondeur d’écoulement est d’au moins 0,1 m. La comparaison est obtenue par l’intermédiaire de l’énoncé conditionnel Con(if_condition, then, else). Les deux rasters (u et h) sont passés comme des objets arcpy.sa.Float() pour s’assurer que le script utilise le type de données pixel correct.
froude = Con(Float(h) > 0.1, Float(u) / SquareRoot(Float(h) * Float(9.81)))Mieux que le critère de 0,1 mètre est de calculer le nombre de Froude partout où la profondeur d’écoulement et la vitesse ont une valeur numérique. arcpy.sa.IsNull() évalue les pixels non numériques. Cependant, nous sommes intéressés par le contraire (c.-à-d., les pixels qui ne sont pas non-numériques), que nous obtenons par le signe ~ (pas) . Le calcul ci-dessous-fonctionnalisé du numéro Froude profite de la possibilité de nid multiple Con() expressions pour vérifier les deux rasters (u et h) pour les pixels numériques.
En outre, nous avons besoin d’une licence d’analyste spatial pour exécuter ce script. Il est donc logique de réécrire le bloc de code ci-dessus dans une fonction qui utilise la fonction spatial_license wrapper de la section checkout Licenses ci-dessus. Pour rendre le code robuste, nous ajoutons également la fonction above-defined *err_info* wrapper.
@err_info
@spatial_license
def calculate_froude(h, u):
return Con(~IsNull(h), Con(~IsNull(u), Float(u) / SquareRoot(Float(h) * Float(9.81))))
froude = calculate_froude(h, u)Pour en savoir plus sur Calculatrice de grille et Map Algebra sur le site Web du développeur (esri)](https://
Statistiques des cellules¶
The evaluation of numerical model data often involves calculating statistical values (e.g., minimum or maximum) of one or more rasters. The comparison of several similar rasters is useful, for example, if the same parameter was calculated with two different models or at different dates (e.g., to assess the morphodynamic evolution of rivers). arcpy.sa.CellStatistics([Raster1, Raster2, ... RasterN], TYPE, ...) enables such statistical evaluations. The following code block illustrates the comparison of flow velocities calculated with two different hydrodynamic numerical models through the calculation of the MEAN (average) and standard deviation (STD).
u_basement = arcpy.Raster("geodata/bm/velocity.tif")
u_tuflow = arcpy.Raster("geodata/tf/velocity.dat")
u_mean = CellStatistics([u_basement, u_tuflow], "MEAN")
u_stdv = CellStatistics([u_basement, u_tuflow], "STD")Lire plus d’options statistiques types et traitement des données non numériques sur le site Web du développeur (Esri).
Opérations de Shapefile¶
Le format vectoriel géospatial shapefile est une invention d’Esri. Pas étonnant, arcpy est bon pour gérer ce format de données vectorielles. Dans le domaine de l’ingénierie hydraulique, cependant, nous créons habituellement des shapefiles (draw) manuellement, soit directement avec ArcGIS ou son concurrent open-source QGIS pour délimiter, par exemple, des régions de débit particulières. Des exemples peuvent être trouvés in the BASEMENT tutorial (explorer la création de points d’élévation, de polygone limite et de fichiers de forme de polyline de rupture). Dans les codes, le traitement des fichiers de forme ne devient important que dans l’analyse de la sortie des modèles numériques (par exemple, pour classer les caractéristiques des unités morphologiques, calculer exactement les zones de patch, ou placer automatiquement les renforts dans les plans de construction). À ce stade, les données raster (sorties de modèles numériques) doivent d’abord être converties en fichiers de forme. C’est pourquoi ce tutoriel commence par la conversion de données de raster en shapefiles ainsi que l’illustration d’autres fonctions telles que le calcul de la surface de patch et l’accès aux tables d’attribut shapefile.
Les rasoirs peuvent être convertis en polygone et en d’autres types de shapefile (par exemple, point). L’exemple suivant présente la conversion d’un raster en un fichier de forme polygone. Il utilise un raster entier de tous les pixels où la profondeur d’écoulement et la vitesse sont inférieures à 1,4 m et 0,15 m/s, respectivement. Ces régions peu profondes et à débit lent sont appelées slackwater (selon Wyrick & Pasternack (2014)). Comme l’habitat de certaines espèces aquatiques est désigné comme étant l’habitat de prédilection de l’eau douce, nous nous demandons maintenant combien de zones d’eau douce le modèle numérique prévoit dans la section simulée de la rivière. À cette fin, nous convertissons le raster en shapefile et calculons la surface du shapefile avec les méthodes arcpy suivantes:
Convert the raster to a shapefile with
arcpy.RasterToPolygon_conversion()with arguments:in_rasterest un objetarcpy.Raster()des valeurs integer (en utilisant unFloatraster entraîne une erreur!).out_polygon_featuresest un string du nom et du répertoire de sortie.simplifyest un string optionnel qui peut être soit"NO_SIMPLIFY"pour forcer le dessin exact des limites du polygone le long de la frontière du pixel, soit"SIMPLIFY"pour permettre les limites du polygone traversant des pixels.
Add a new field to the new polygon shapefile with
arcpy.AddField_management()with arguments:field_namepeut être n’importe quelle chaîne sans blanc.field_typeest un string définissant si le champ est numérique (par exemple,"FLOAT"ou"LONG"pour integer), date/heure ("DATE"),"TEXT", ou"RASTER".field_precisionest un (facultatif) integer (Long) définissant le nombre de chiffres qui peuvent être stockés dans le nouveau champ.D’autres arguments optionnels peuvent être définis pour définir le nombre de décimales ou de caractères, un nom de champ alternatif, activer
NULL, un domaine de champ, ou si un champ est requis.
Calculate patch area with
arcpy.CalculateGeometryAttributes_management()with arguments:in_featuresest un string définissant le répertoire et le nom d’un calque de fonction.geometry_propertyest une liste imbriquée de[[Target-field-name, Property], [Another-target-field-name, Another-property], ...]pour calculer des propriétés géométriques telles que"AREA","HOLE_COUNT", ou"PART_COUNT"(et beaucoup d’autres).area_unitpeut être"SQUARE_METERS"ou"SQUARE_KILOMETERS"(et beaucoup d’autres options pour les unités habituelles des États-Unis).D’autres arguments optionnels peuvent être définis pour définir des unités de longueur (pour les évaluations du périmètre), ou un système de coordonnées et un format.
Le bloc de code ci-dessous présente en outre l’application de ces méthodes et montre également comment la valeur de zone peut être lue ligne par ligne à partir de la table d’attribut d’un fichier de forme avec arcpy.da.UpdateCursor(shapefile-name, field-name).
# create a slackwater raster that is arcpy.sa.Int(1) where h and u criteria are true and NULL elsewhere
slackwater = Con((Float(h) <= 1.4) & (Float(u) <= 0.15), Int(1))
# define directory and name of the new shapefile
new_shp_file = "geodata/shapefiles/slackwater.shp"
# convert slackwater raster to polygon shapefile
arcpy.RasterToPolygon_conversion(in_raster=slackwater, out_polygon_features=new_shp_file, simplify="NO_SIMPLIFY")
# add a new field to the new shapefile's attribute table (name)
arcpy.AddField_management(new_shp_file, field_name="F_AREA", field_type="FLOAT", field_precision=9)
# calculate area of all polygons in attribute table
area_unit="SQUARE_METERS"
arcpy.CalculateGeometryAttributes_management(in_features=new_shp_file, geometry_property=[["F_AREA", "AREA"]], area_unit=area_unit)
area = 0.0
with arcpy.da.UpdateCursor(new_shp_file, "F_AREA") as cursor:
for row in cursor:
try:
area += float(row[0])
except ValueError:
print("WARNING: Patch with invalid area value (%s)." % str(row))
print("Sum of all patches = {0} {1}".format(str(area), area_unit))- Wyrick, J. R., & Pasternack, G. B. (2014). Geospatial organization of fluvial landforms in a gravel-cobble river: Beyond the riffle-pool couplet. Geomorphology, 213(Supplement C), 48–65. 10.1016/j.geomorph.2013.12.040