Est-il possible de regarder le contenu de Shapefile en utilisant Python sans licence ArcMap?

40

Je me suis demandé s'il était possible d'examiner le contenu d'un fichier de formes à l'aide de Python sans licence ArcMap. La situation est que vous pouvez créer des fichiers de formes à partir de nombreuses applications différentes, pas seulement à partir du logiciel ESRI. J'aimerais créer un script Python qui vérifie la référence spatiale, le type d'entité, les noms d'attributs et les définitions, ainsi que le contenu des champs d'un fichier de formes, et les compare à un ensemble de valeurs acceptables. J'aimerais que ce script fonctionne même si l'organisation n'a pas de licence ESRI. Pour faire quelque chose comme ceci, devez-vous utiliser ArcPy ou pouvez-vous creuser un fichier de formes sans utiliser ArcPy?

python shapefile Caitlin
la source
1
Cela dépend de l'effort que vous voulez y mettre. Il y a plusieurs bibliothèques open source qui aideront (j'aime OGR selon la réponse d'Aarons), mais si vous voulez vraiment contrôler (et êtes prêt à y travailler), le Shapefile (à l'origine par Esri) est un format ouvert, voir en.wikipedia.org/wiki/Shapefile
Michael Stimson
1
Les fichiers de formes ESRI récents (les deux dernières années) sont cachés dans leur nouveau format de géodatabase. Il semble que rien ne puisse les casser sauf le logiciel ARCxxx. Beaucoup d'organismes publics l'utilisent pour informer le public ... dommage.

Réponses:

34

Je vous recommanderais de vous familiariser avec l' API Python GDAL / OGR pour pouvoir utiliser à la fois des données vectorielles et raster. Le moyen le plus simple de commencer à utiliser GDAL / OGR consiste à utiliser une distribution python telle que python (x, y) , Anaconda ou OSGeo4W .

Plus de détails sur l'utilisation de GDAL pour vos tâches spécifiques:

De plus, je recommanderais le tutoriel suivant de USU pour vous aider à démarrer.

Reprenant les exemples ci-dessus, le script suivant utilise les outils FOSS pour effectuer les opérations suivantes:

  1. Vérifier la référence spatiale
  2. Obtenir des champs et des types de fichiers de formes
  3. Vérifier si les lignes d'un champ défini par l'utilisateur contiennent une valeur
# Import the necessary modules
from  osgeo import ogr, osr

driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open(r'C:\your\shapefile.shp')

# Get Projection from layer
layer = shp.GetLayer()
spatialRef = layer.GetSpatialRef()
print spatialRef

# Get Shapefile Fields and Types
layerDefinition = layer.GetLayerDefn()

print "Name  -  Type  Width  Precision"
for i in range(layerDefinition.GetFieldCount()):
    fieldName =  layerDefinition.GetFieldDefn(i).GetName()
    fieldTypeCode = layerDefinition.GetFieldDefn(i).GetType()
    fieldType = layerDefinition.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
    fieldWidth = layerDefinition.GetFieldDefn(i).GetWidth()
    GetPrecision = layerDefinition.GetFieldDefn(i).GetPrecision()
    print fieldName + " - " + fieldType+ " " + str(fieldWidth) + " " + str(GetPrecision)

# Check if rows in attribute table meet some condition
inFeature = layer.GetNextFeature()
while inFeature:

    # get the cover attribute for the input feature
    cover = inFeature.GetField('cover')

    # check to see if cover == grass
    if cover == 'trees':
        print "Do some action..."

    # destroy the input feature and get a new one
    inFeature = None
    inFeature = inLayer.GetNextFeature()
Aaron
la source
Bien que ce ne soit pas un gros problème, certains objets contiennent une méthode Destroy (), mais vous ne devriez jamais l'utiliser
Mike T
Merci pour la perspicacité @MikeT. La documentation GDAL / OGR utilise la méthode 'Destroy ()' tout au long de leur livre de recettes. Quels problèmes voyez-vous avec cette méthode?
Aaron
1
Il existe des situations où des erreurs de segmentation peuvent se produire lorsque vous utilisez Destroy (), et exposer cette méthode dans les liaisons était une erreur de conception. Une meilleure méthode consiste à déréférencer des objets GDAL comme inFeature = None. Le livre de recettes GDAL / OGR ne fait pas partie de, ni n'a été écrit par l'équipe principale de GDAL / OGR.
Mike T
@MikeT J'ai modifié le message pour inclure vos commentaires - merci.
Aaron
31

Il existe de nombreux modules pour lire les fichiers de formes en Python, plus anciens que ArcPy, consultez l’ index des packages Python (PyPi): fichiers de formes . Il y a aussi beaucoup d'exemples dans GIS SE (recherche de [Python] Fiona , par exemple)

Tous peuvent lire la géométrie, les champs et les projections.

Mais d’autres modules comme PySAL: la bibliothèque d’analyses spatiales Python , Cartopy (qui utilisent pyshp ) ou Matplotlib Basemap peuvent également lire des fichiers de formes, entre autres.

Le plus simple à utiliser est Fiona , mais si vous ne connaissez que ArcPy, utilisez pyshp , car osgeo et Fiona exigent l’ installation de la bibliothèque GDAL C / C ++, GeoPandas nécessite le module Pandas et PySAL est trop volumineux (beaucoup d’autres traitements).

Si vous souhaitez uniquement lire le contenu d'un fichier de formes, vous n'avez pas besoin de choses complexes, utilisez simplement le protocole d' interface géographique (GeoJSON) également implémenté dans ArcPy ( ArcPy: AsShape ).

Avec Fiona (en tant que dictionnaires Python):

import fiona
with fiona.open('a_shape.shp') as shp:
     # schema of the shapefile
     print shp.schema
     {'geometry': 'Point', 'properties': OrderedDict([(u'DIP', 'int:2'), (u'DIP_DIR', 'int:3'), (u'TYPE', 'str:10')])}
     # projection
     print shp.crs
     {u'lon_0': 4.367486666666666, u'ellps': u'intl', u'y_0': 5400088.438, u'no_defs': True, u'proj': u'lcc', u'x_0': 150000.013, u'units': u'm', u'lat_2': 49.8333339, u'lat_1': 51.16666723333333, u'lat_0': 90}
     for feature in shp:
        print feature              
{'geometry': {'type': 'Point', 'coordinates': (272070.600041, 155389.38792)}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'DIP', 30), (u'DIP_DIR', 130), (u'TYPE', u'incl')])}
{'geometry': {'type': 'Point', 'coordinates': (271066.032148, 154475.631377)}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'DIP', 55), (u'DIP_DIR', 145), (u'TYPE', u'incl')])}
{'geometry': {'type': 'Point', 'coordinates': (273481.498868, 153923.492988)}, 'type': 'Feature', 'id': '2', 'properties': OrderedDict([(u'DIP', 40), (u'DIP_DIR', 155), (u'TYPE', u'incl')])}

Avec pyshp (sous forme de dictionnaires Python)

import shapefile
reader= shapefile.Reader("a_shape.shp")
# schema of the shapefile
print dict((d[0],d[1:]) for d in reader.fields[1:])
{'DIP_DIR': ['N', 3, 0], 'DIP': ['N', 2, 0], 'TYPE': ['C', 10, 0]}
fields = [field[0] for field in reader.fields[1:]]
for feature in reader.shapeRecords():
    geom = feature.shape.__geo_interface__
    atr = dict(zip(fields, feature.record))
    print geom, atr
{'type': 'Point', 'coordinates': (272070.600041, 155389.38792)} {'DIP_DIR': 130, 'DIP': 30, 'TYPE': 'incl'}
{'type': 'Point', 'coordinates': (271066.032148, 154475.631377)} {'DIP_DIR': 145, 'DIP': 55, 'TYPE': 'incl'}
{'type': 'Point', 'coordinates': (273481.498868, 153923.492988)} {'DIP_DIR': 155, 'DIP': 40, 'TYPE': 'incl'}

Avec osgeo / ogr (sous forme de dictionnaires Python)

from osgeo import ogr
reader = ogr.Open("a_shape.shp")
layer = reader.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    print feature.ExportToJson()
{"geometry": {"type": "Point", "coordinates": [272070.60004, 155389.38792]}, "type": "Feature", "properties": {"DIP_DIR": 130, "DIP": 30, "TYPE": "incl"}, "id": 0}
{"geometry": {"type": "Point", "coordinates": [271066.032148, 154475.631377]}, "type": "Feature", "properties": {"DIP_DIR": 145, "DIP": 55, "TYPE": "incl"}, "id": 1}
{"geometry": {"type": "Point", "coordinates": [273481.49887, 153923.492988]}, "type": "Feature", "properties": {"DIP_DIR": 155, "DIP": 40, "TYPE": "incl"}, "id": 2}

Avec GeoPandas (en tant que cadre de données Pandas)

import geopandas as gp
shp = gp.GeoDataFrame.from_file('a_shape.shp')
print shp
        DIP_DIR    DIP  TYPE                       geometry
0         130       30  incl          POINT (272070.600041 155389.38792)
1         145       55  incl          POINT (271066.032148 154475.631377)
2         155       40  incl          POINT (273481.498868 153923.492988)

* note sur les géopandas Vous devez utiliser les anciennes versions de Fiona et GDAL avec celle-ci, sinon l'installation ne s'installera pas. GDAL: 1.11.2 Fiona: 1.6.0 Géopandas: 0.1.0.dev-

Il existe de nombreux tutoriels sur le Web et même des livres ( développement géospatial Python , apprentissage de l'analyse géospatiale avec Python et géotraitement avec Python , sous presse).

Plus généralement, si vous souhaitez utiliser Python sans ArcPy, reportez-vous à la rubrique Mappage thématique simple du fichier de formes à l'aide de Python?

gène
la source
Notez que la page principale de Fiona indiqueThe kinds of data in GIS are roughly divided into rasters representing continuous scalar fields (land surface temperature or elevation, for example) and vectors representing discrete entities like roads and administrative boundaries. Fiona is concerned exclusively with the latter
Mawg le
2
Évidemment, la question concerne les fichiers de formes et non les rasters. Ce sont d'autres modules pour les fichiers raster.
gène le
Très bonne réponse! Quelque chose à mettre à jour en 2017?
Michael