Is it possible to look at the contents of Shapefile using Python without an ArcMap license?

Question

I wondered if it is possible to look at the contents of a shapefile using Python without having and ArcMap license. The situation is that you can create shapefiles from many different applications, not only from ESRI software. I would like to create a Python script that checks the spatial reference, feature type, attribute names and definitions, and the contents of the fields in a shapefile and compares them to a set of acceptable values. I would like this script to work even if the organization does not have any ESRI licenses. To do something like this, do you have to use ArcPy or can you dig into a shapefile without using ArcPy?

Answer 1

I would recommend becoming familiar with the Python GDAL/OGR API to work with both vector and raster data. The easiest way to start using GDAL/OGR is via a python distribution such as python(x,y), Anaconda, or OSGeo4W.

Further details on using GDAL for your specific tasks:

• Get Shapefile Fields and Types
• Get Projection

Additionally, I would recommend the following tutorial from USU to get you started.

• Geoprocessing with Python using Open Source GIS

---

Borrowing from the examples above, the following script uses FOSS tools to perform the following actions:

1. Check the spatial reference
2. Get shapefile fields and types
3. Check if rows in a user-defined field contain some value

---

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

Answer 2

There are many modules to read shapefiles in Python, older than ArcPy, look at the Python Package Index (PyPi): shapefiles. There are also many examples in GIS SE (search for [Python] Fiona, for example)

All can read the geometry, the fields and the projections.

• The older is osgeo (GDAL/OGR), look at the Python GDAL/OGR Cookbook for example
• another solution is Fiona, also based on GDAL/OGR, but easier to use (with Python dictionaries: GeoJSON format).
• pyshp (shapefile) is a pure Python solution
• GeoPandas use Fiona to read/write the shapefiles and Pandas, for data analysis tools. You need to know Pandas to use it.

But other modules as PySAL:the Python Spatial Analysis Library, Cartopy (which use pyshp) or Matplotlib Basemap can also read shapefiles, among other things.

The easiest to use is Fiona, but if you only know ArcPy, use pyshp, because osgeo and Fiona require that the GDAL C/C++ library be installed, GeoPandas needs the Pandas module and PySAL is too big (many, many others treatments)

If you only want to read the content of a shapefile, you don't need complex things, simply use the geo interface protocol (GeoJSON) also implemented in ArcPy (ArcPy: AsShape)

With Fiona (as Python dictionaries):

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

With pyshp (as Python dictionaries)

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

With osgeo/ogr (as Python dictionaries)

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}

With GeoPandas (as Pandas dataframe)

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 on geopandas You have to use older versions of Fiona and GDAL with it or it won't install. GDAL: 1.11.2 Fiona: 1.6.0 Geopandas: 0.1.0.dev-

There are many tutorials on the Web and even books (Python Geospatial Development , Learning Geospatial Analysis with Python and Geoprocessing with Python, in press)

More generally, if you want to use Python without ArcPy, look at Simple thematic mapping of shapefile using Python?

Answer 3

There are geospatial Python libraries besides ArcPy that will give you these abilities. Here are two examples:

The Python Shapefile Library (pyshp)

GeoPandas

If you're interested in other libraries, this post on essential Python Geospatial libraries is a good place to look.