Ist es möglich, den Inhalt von Shapefile mit Python ohne eine ArcMap-Lizenz anzuzeigen?

40

Ich habe mich gefragt, ob es möglich ist, den Inhalt eines Shapefiles mit Python zu betrachten, ohne eine ArcMap-Lizenz zu besitzen. Die Situation ist, dass Sie Shapefiles aus vielen verschiedenen Anwendungen erstellen können, nicht nur aus ESRI-Software. Ich möchte ein Python-Skript erstellen, das den Raumbezug, den Feature-Typ, die Attributnamen und -definitionen sowie den Inhalt der Felder in einem Shapefile überprüft und mit einer Reihe akzeptabler Werte vergleicht. Ich möchte, dass dieses Skript auch dann funktioniert, wenn die Organisation keine ESRI-Lizenzen hat. Müssen Sie dafür ArcPy verwenden oder können Sie in ein Shapefile graben, ohne ArcPy zu verwenden?

Caitlin
quelle
1
Es hängt davon ab, wie viel Aufwand Sie betreiben möchten. Es gibt mehrere Open-Source-Bibliotheken, die helfen (ich mag OGR gemäß Aarons Antwort), aber wenn Sie wirklich die Kontrolle über das Shapefile haben möchten (und bereit sind, dafür zu arbeiten) (ursprünglich von Esri) ist ein offenes Format, siehe en.wikipedia.org/wiki/Shapefile
Michael Stimson
1
Aktuelle (letzte paar Jahre) ESRI-Shapefiles sind in ihrem neuen Geodatabase-Format versteckt. Es scheint, dass nichts sie außer ARCxxx-Software brechen kann. Viele öffentliche Stellen nutzen es zur Information der Öffentlichkeit ... schade.

Antworten:

34

Ich würde empfehlen, sich mit der Python GDAL / OGR- API vertraut zu machen , um sowohl mit Vektor- als auch mit Rasterdaten arbeiten zu können. Der einfachste Weg, um mit GDAL / OGR zu beginnen, ist über eine Python-Distribution wie Python (x, y) , Anaconda oder OSGeo4W .

Weitere Details zur Verwendung von GDAL für Ihre spezifischen Aufgaben:

Zusätzlich würde ich das folgende Tutorial von USU empfehlen, um Ihnen den Einstieg zu erleichtern.


In Anlehnung an die obigen Beispiele verwendet das folgende Skript FOSS-Tools, um die folgenden Aktionen auszuführen:

  1. Überprüfen Sie den Raumbezug
  2. Ruft Shapefile-Felder und -Typen ab
  3. Überprüfen Sie, ob Zeilen in einem benutzerdefinierten Feld einen Wert enthalten

# 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
quelle
Danke für den Einblick @MikeT. In der GDAL / OGR-Dokumentation wird im gesamten Kochbuch die Methode 'Destroy ()' verwendet. Welche Probleme sehen Sie mit dieser Methode?
Aaron
1
Es gibt Situationen, in denen Segfaults auftreten können, wenn Sie Destroy () verwenden, und es war ein Entwurfsfehler, diese Methode in den Bindungen verfügbar zu machen. Eine bessere Methode ist die Dereferenzierung von GDAL-Objekten wie inFeature = None. Das GDAL / OGR-Kochbuch ist nicht Teil des GDAL / OGR-Kernteams oder wird von diesem geschrieben.
Mike T
@MikeT Ich habe den Beitrag so bearbeitet, dass er Ihre Kommentare enthält - danke.
Aaron
31

Es gibt viele Module zum Lesen von Shapefiles in Python, die älter als ArcPy sind. Schauen Sie sich den Python Package Index (PyPi) an: Shapefiles . Es gibt auch viele Beispiele in GIS SE (suchen Sie beispielsweise nach [Python] Fiona )

Alle können die Geometrie, die Felder und die Projektionen lesen.

Aber auch andere Module wie PySAL: Die Python Spatial Analysis Library , Cartopy (die pyshp verwenden ) oder Matplotlib Basemap können unter anderem Shapefiles lesen.

Am einfachsten zu bedienen ist Fiona , aber wenn Sie nur ArcPy kennen, verwenden Sie pyshp , da osgeo und Fiona die Installation der GDAL C / C ++ - Bibliothek erfordern , GeoPandas das Pandas- Modul benötigt und PySAL zu groß ist (viele, viele andere Behandlungen).

Wenn Sie nur den Inhalt eines Shapefiles lesen möchten, brauchen Sie keine komplexen Dinge. Verwenden Sie einfach das Geo-Interface- Protokoll (GeoJSON), das auch in ArcPy ( ArcPy: AsShape ) implementiert ist.

Mit Fiona (als Python-Wörterbücher):

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')])}

Mit pyshp (als Python-Wörterbücher)

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'}

Mit osgeo / ogr (als Python-Wörterbücher)

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}

Mit GeoPandas (als Pandas-Datenrahmen)

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)

* Hinweis zu Geopandas Sie müssen ältere Versionen von Fiona und GDAL verwenden, sonst wird es nicht installiert. GDAL: 1.11.2 Fiona: 1.6.0 Geopandas: 0.1.0.dev-

Es gibt viele Tutorials im Web und sogar Bücher ( Python Geospatial Development , Learning Geospatial Analysis mit Python und Geoprocessing mit Python , im Druck)

Wenn Sie Python ohne ArcPy verwenden möchten, lesen Sie allgemein den Abschnitt Einfache thematische Zuordnung von Shapefile mit Python.

Gen
quelle
Beachten Sie, dass auf der Hauptseite von Fiona stehtThe 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
2
Offensichtlich geht es um Shapefiles und nicht um Raster. Sie sind andere Module für Rasterdateien.
Gen
Gute Antwort! Gibt es 2017 etwas zu aktualisieren?
Michael