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?
-
1It 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
-
1Recent (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
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:
- Check the spatial reference
- Get shapefile fields and types
- 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()
-
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
-
1There 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 -
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?
-
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 says reinstate Monica Sep 21 '16 at 9:52 -
2Evident, the question is about shapefiles and not rasters. They are other modules for raster files. – gene Sep 21 '16 at 15:19
-
There are geospatial Python libraries besides ArcPy that will give you these abilities. Here are two examples:
The Python Shapefile Library (pyshp)
If you're interested in other libraries, this post on essential Python Geospatial libraries is a good place to look.