Czy można przeglądać zawartość Shapefile przy użyciu Pythona bez licencji ArcMap?

40

Zastanawiałem się, czy można przeglądać zawartość pliku shapefile przy użyciu Pythona bez posiadania licencji ArcMap. Sytuacja polega na tym, że możesz tworzyć pliki shapefile z wielu różnych aplikacji, nie tylko z oprogramowania ESRI. Chciałbym utworzyć skrypt w języku Python, który sprawdza odniesienie przestrzenne, typ funkcji, nazwy atrybutów i definicje oraz zawartość pól w pliku kształtu i porównuje je z zestawem dopuszczalnych wartości. Chciałbym, aby ten skrypt działał, nawet jeśli organizacja nie ma żadnych licencji ESRI. Aby zrobić coś takiego, czy musisz użyć ArcPy, czy możesz zagłębić się w plik kształtu bez korzystania z ArcPy?

Caitlin
źródło
1
To zależy od tego, ile wysiłku chcesz w to włożyć .. istnieje kilka bibliotek open source, które pomogą (lubię OGR zgodnie z odpowiedzią Aaronsa), ale jeśli naprawdę chcesz kontrolować (i jesteś gotowy na to pracować), Shapefile (pierwotnie Esri) jest formatem otwartym, patrz en.wikipedia.org/wiki/Shapefile
Michael Stimson
1
Ostatnie (ostatnie kilka lat) pliki kształtów ESRI są ukryte w nowym formacie geobazy. Wygląda na to, że nic nie może ich złamać oprócz oprogramowania ARCxxx. Wiele agencji publicznych używa go do informacji publicznej ... wstyd.

Odpowiedzi:

34

Polecam zapoznanie się z interfejsem API Python GDAL / OGR do pracy z danymi wektorowymi i rastrowymi. Najłatwiejszym sposobem na rozpoczęcie korzystania z GDAL / OGR jest dystrybucja pytona, taka jak python (x, y) , Anaconda lub OSGeo4W .

Dalsze szczegóły dotyczące korzystania z GDAL do określonych zadań:

Ponadto polecam następujący samouczek od USU, aby zacząć.


Pożyczając z powyższych przykładów, poniższy skrypt używa narzędzi FOSS do wykonywania następujących czynności:

  1. Sprawdź odniesienie przestrzenne
  2. Uzyskaj pola i typy plików shapefile
  3. Sprawdź, czy wiersze w polu zdefiniowanym przez użytkownika zawierają jakąś wartość

# 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
źródło
Dzięki za wgląd @MikeT. Dokumentacja GDAL / OGR używa metody „Destroy ()” w całej książce kucharskiej. Jakie problemy widzisz przy tej metodzie?
Aaron
1
Istnieją sytuacje, w których mogą wystąpić awarie segfault podczas korzystania z funkcji Destroy (), a ujawnienie tej metody w powiązaniach było błędem projektowym. Lepszą metodą jest usunięcie dereferencji z obiektów GDAL, takich jak inFeature = None. Książka kucharska GDAL / OGR nie jest częścią ani nie jest napisana przez podstawowy zespół GDAL / OGR.
Mike T
@MikeT Zredagowałem ten post, aby dodać komentarze - dzięki.
Aaron
31

Istnieje wiele modułów do odczytu plików kształtów w Pythonie, starszych niż ArcPy, spójrz na Indeks pakietów Python (PyPi): pliki kształtów . Istnieje również wiele przykładów w GIS SE (na przykład wyszukaj [Python] Fiona )

Wszyscy mogą odczytać geometrię, pola i rzuty.

Ale inne moduły, takie jak PySAL: biblioteka Python Spatial Analysis Library , Cartopy (które używają pyshp ) lub Matplotlib Basemap mogą również odczytywać pliki shapefile .

Najłatwiejszy w użyciu jest Fiona , ale jeśli znasz tylko ArcPy, użyj pyshp , ponieważ osgeo i Fiona wymagają zainstalowania biblioteki GDAL C / C ++, GeoPandas potrzebuje modułu Pandas , a PySAL jest zbyt duży (wiele, wiele innych zabiegów)

Jeśli chcesz tylko czytać zawartość pliku shapefile, nie potrzebujesz skomplikowanych rzeczy, po prostu użyj protokołu interfejsu geo (GeoJSON) również zaimplementowanego w ArcPy ( ArcPy: AsShape )

Z Fiona (jako słowniki 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')])}

Z pyshp (jako słowniki 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'}

Z osgeo / ogr (jako słowniki 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}

Z GeoPandas (jako ramka danych 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)

* uwaga na geopandas Musisz używać z nim starszych wersji Fiony i GDAL, inaczej się nie zainstaluje. GDAL: 1.11.2 Fiona: 1.6.0 Geopandy: 0.1.0.dev-

W sieci istnieje wiele samouczków, a nawet książek ( Python Geospatial Development , Learning Geospatial Analysis with Python i Geoprocessing with Python , in press)

Mówiąc bardziej ogólnie, jeśli chcesz używać Pythona bez ArcPy, spójrz na Proste tematyczne mapowanie pliku kształtu za pomocą Pythona?

gen
źródło
Zauważ, że strona główna Fiony mówiThe 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
Oczywiste jest, że pytanie dotyczy plików kształtowych, a nie rastrów. Są to inne moduły dla plików rastrowych.
gen
Świetna odpowiedź! Coś do zaktualizowania w 2017 roku?
Michael