Apakah mungkin untuk melihat konten Shapefile menggunakan Python tanpa lisensi ArcMap?

Saya bertanya-tanya apakah mungkin untuk melihat isi dari suatu shapefile menggunakan Python tanpa memiliki dan lisensi ArcMap. Situasinya adalah Anda dapat membuat shapefile dari berbagai aplikasi, tidak hanya dari perangkat lunak ESRI. Saya ingin membuat skrip Python yang memeriksa referensi spasial, tipe fitur, nama atribut dan definisi, dan isi bidang dalam sebuah shapefile dan membandingkannya dengan seperangkat nilai yang dapat diterima. Saya ingin agar skrip ini berfungsi bahkan jika organisasi tidak memiliki lisensi ESRI. Untuk melakukan sesuatu seperti ini, apakah Anda harus menggunakan ArcPy atau dapatkah Anda menggali ke dalam shapefile tanpa menggunakan ArcPy?

Saya akan merekomendasikan untuk menjadi terbiasa dengan API Python GDAL / OGR untuk bekerja dengan data vektor dan raster. Cara termudah untuk mulai menggunakan GDAL / OGR adalah melalui distribusi python seperti python (x, y) , Anaconda , atau OSGeo4W .

Rincian lebih lanjut tentang penggunaan GDAL untuk tugas spesifik Anda:

• Dapatkan Bidang dan Jenis Shapefile
• Dapatkan Proyeksi

Selain itu, saya akan merekomendasikan tutorial berikut dari USU untuk membantu Anda memulai.

• Geoprocessing dengan Python menggunakan Open Source GIS

Meminjam dari contoh di atas, skrip berikut menggunakan alat FOSS untuk melakukan tindakan berikut:

1. Periksa referensi spasial
2. Dapatkan bidang dan jenis shapefile
3. Periksa apakah baris dalam bidang yang ditentukan pengguna berisi beberapa nilai

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

Ada banyak modul untuk membaca shapefile dengan Python, lebih tua dari ArcPy, lihat pada Python Package Index (PyPi): shapefile . Ada juga banyak contoh di GIS SE (mencari [Python] Fiona , misalnya)

Semua dapat membaca geometri, bidang, dan proyeksi.

• Yang lebih tua adalah osgeo (GDAL / OGR) , lihat Cookbook Python GDAL / OGR misalnya
• solusi lain adalah Fiona , juga didasarkan pada GDAL / OGR, tetapi lebih mudah digunakan (dengan kamus Python: format GeoJSON).
• pyshp (shapefile) adalah solusi Python murni
• GeoPandas menggunakan Fiona untuk membaca / menulis shapefile dan Pandas , untuk alat analisis data. Anda perlu tahu Panda untuk menggunakannya.

Tetapi modul-modul lain seperti PySAL: Python Spatial Analysis Library , Cartopy (yang menggunakan pyshp ) atau Matplotlib Basemap juga dapat membaca shapefile, antara lain.

Cara termudah untuk digunakan adalah Fiona , tetapi jika Anda hanya tahu ArcPy, gunakan pyshp , karena osgeo dan Fiona mengharuskan pustaka GDAL C / C ++ diinstal, GeoPandas membutuhkan modul Pandas dan PySAL terlalu besar (banyak, banyak perawatan lain)

Jika Anda hanya ingin membaca konten shapefile, Anda tidak perlu hal-hal yang kompleks, cukup gunakan protokol antarmuka geo (GeoJSON) juga diimplementasikan dalam ArcPy ( ArcPy: AsShape )

Dengan Fiona (sebagai kamus 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')])}

Dengan pyshp (seperti kamus 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'}

Dengan osgeo / ogr (sebagai kamus 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}

Dengan GeoPandas (sebagai kerangka data 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)

* catatan di geopanda Anda harus menggunakan Fiona dan GDAL versi lama dengan itu atau itu tidak akan menginstal. GDAL: 1.11.2 Fiona: 1.6.0 Geopanda: 0.1.0.dev-

Ada banyak tutorial di Web dan bahkan buku-buku ( Pengembangan Geospasial Python , Belajar Analisis Geospasial dengan Python dan Geoprocessing dengan Python , sedang dicetak )

Lebih umum, jika Anda ingin menggunakan Python tanpa ArcPy, lihat pemetaan tematik sederhana dari shapefile menggunakan Python?

Ada perpustakaan Python geospasial selain ArcPy yang akan memberi Anda kemampuan ini. Berikut adalah dua contoh:

Perpustakaan Python Shapefile (pyshp)

GeoPandas

Jika Anda tertarik pada pustaka lain, posting tentang pustaka Python Geospatial yang esensial ini adalah tempat yang baik untuk melihat.