40

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?

|improve this question
  • 1
    It depends on how much effort you want to put into it.. there are several open source libraries that will help (I like OGR as per Aarons' answer) but if you really want control (and are prepared to work for it) the Shapefile (originally by Esri) is an open format see en.wikipedia.org/wiki/Shapefile – Michael Stimson May 5 '15 at 2:27
  • 1
    Recent (last couple of years) ESRI shapefiles are hidden in their new geodatabase format. It seems that nothing can break them except ARCxxx software. Many public agencies are using it for public information...shame. – user51663 May 6 '15 at 15:45
34

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:

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

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()
|improve this answer
  • While not a big issue, certain objects contain a Destroy() method, but you should never use it – Mike T Dec 13 '16 at 20:24
  • Thanks for the insight @MikeT. The GDAL/OGR documentation uses the 'Destroy()' method throughout their cookbook. What issues do you see with that method? – Aaron♦ Dec 14 '16 at 17:28
  • 1
    There are situations where segfaults may happen when you use Destroy(), and it was a design mistake to expose this method in the bindings. A better method is to dereference GDAL objects like inFeature = None. The GDAL/OGR cookbook is not part of, or written by the GDAL/OGR core team. – Mike T Dec 14 '16 at 23:12
  • @MikeT I've edited the post to include your comments--thanks. – Aaron♦ Dec 15 '16 at 20:04
31

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.

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?

|improve this answer
  • Note that Fiona's main page says The 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 Sep 21 '16 at 9:52
  • 2
    Evident, the question is about shapefiles and not rasters. They are other modules for raster files. – gene Sep 21 '16 at 15:19
  • Great answer! Anything to update in 2017? – Michael Dec 2 '17 at 7:56
11

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.

|improve this answer
  • 1
    +1 for pyshp, super easy and straightforward to use. Sounds like everything the OP wants would be accommodated. – mr.adam May 4 '15 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.