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

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?

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ń:

• Uzyskaj pola i typy Shapefile
• Uzyskaj projekcję

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

• Geoprocessing w Pythonie przy użyciu GIS typu Open Source

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()

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.

• Starszym jest osgeo (GDAL / OGR) , patrz na przykład książka kucharska Python GDAL / OGR
• innym rozwiązaniem jest Fiona , również oparta na GDAL / OGR, ale łatwiejsza w użyciu (ze słownikami Python: format GeoJSON).
• pyshp (shapefile) to czyste rozwiązanie Pythona
• GeoPandy używają Fiony do odczytu / zapisu plików kształtów i Pand dla narzędzi analizy danych. Musisz znać Pandy, aby z niego korzystać.

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?

Oprócz ArcPy istnieją geoprzestrzenne biblioteki Pythona, które dają te możliwości. Oto dwa przykłady:

Biblioteka kształtów Python (pyshp)

GeoPandy

Jeśli interesują Cię inne biblioteki, ten post na temat podstawowych bibliotek geoprzestrzennych Pythona jest dobrym miejscem do obejrzenia.