The Python Shapefile Library (pyshp) reads and writes ESRI Shapefiles in pure Python.
The Python Shapefile Library (pyshp) provides read and write support for the Esri Shapefile format. The Shapefile format is a popular Geographic Information System vector data format created by Esri. For more information about this format please read the well-written "ESRI Shapefile Technical Description - July 1998" located at http://www.esri.com/library/whitepapers/p dfs/shapefile.pdf . The Esri document describes the shp and shx file formats. However a third file format called dbf is also required. This format is documented on the web as the "XBase File Format Description" and is a simple file-based database format created in the 1960's. For more on this specification see: [http://www.clicketyclick.dk/databases/xbase/format/index.html](http://www.clicketyclick.d k/databases/xbase/format/index.html)
Both the Esri and XBase file-formats are very simple in design and memory efficient which is part of the reason the shapefile format remains popular despite the numerous ways to store and exchange GIS data available today.
Pyshp is compatible with Python 2.7-3.x.
This document provides examples for using pyshp to read and write shapefiles. However many more examples are continually added to the pyshp wiki on GitHub, the blog http://GeospatialPython.com, and by searching for pyshp on https://gis.stackexchange.com.
Currently the sample census blockgroup shapefile referenced in the examples is available on the GitHub project site at https://github.com/GeospatialPython/pyshp. These examples are straight-forward and you can also easily run them against your own shapefiles with minimal modification.
Important: If you are new to GIS you should read about map projections. Please visit: https://github.com/GeospatialPython/pyshp/wiki/Map-Projections
I sincerely hope this library eliminates the mundane distraction of simply reading and writing data, and allows you to focus on the challenging and FUN part of your geospatial project.
Before doing anything you must import the library.
>>> import shapefile
The examples below will use a shapefile created from the U.S. Census Bureau Blockgroups data set near San Francisco, CA and available in the git repository of the pyshp GitHub site.
To read a shapefile create a new "Reader" object and pass it the name of an existing shapefile. The shapefile format is actually a collection of three files. You specify the base filename of the shapefile or the complete filename of any of the shapefile component files.
>>> sf = shapefile.Reader("shapefiles/blockgroups")
OR
>>> sf = shapefile.Reader("shapefiles/blockgroups.shp")
OR
>>> sf = shapefile.Reader("shapefiles/blockgroups.dbf")
OR any of the other 5+ formats which are potentially part of a shapefile. The library does not care about file extensions.
You can also load shapefiles from any Python file-like object using keyword arguments to specify any of the three files. This feature is very powerful and allows you to load shapefiles from a url, from a zip file, serialized object, or in some cases a database.
>>> myshp = open("shapefiles/blockgroups.shp", "rb")
>>> mydbf = open("shapefiles/blockgroups.dbf", "rb")
>>> r = shapefile.Reader(shp=myshp, dbf=mydbf)
Notice in the examples above the shx file is never used. The shx file is a very simple fixed-record index for the variable length records in the shp file. This file is optional for reading. If it's available pyshp will use the shx file to access shape records a little faster but will do just fine without it.
Shapefiles have a number of attributes for inspecting the file contents. A shapefile is a container for a specific type of geometry, and this can be checked using the shapeType attribute.
>>> sf.shapeType
5
Shape types are represented by numbers between 0 and 31 as defined by the shapefile specification and listed below. It is important to note that numbering system has several reserved numbers which have not been used yet therefore the numbers of the existing shape types are not sequential:
- NULL = 0
- POINT = 1
- POLYLINE = 3
- POLYGON = 5
- MULTIPOINT = 8
- POINTZ = 11
- POLYLINEZ = 13
- POLYGONZ = 15
- MULTIPOINTZ = 18
- POINTM = 21
- POLYLINEM = 23
- POLYGONM = 25
- MULTIPOINTM = 28
- MULTIPATCH = 31
Based on this we can see that our blockgroups shapefile contains Polygon type shapes. The shape types are also defined as constants in the shapefile module, so that we can compare types more intuitively:
>>> sf.shapeType == shapefile.POLYGON
True
Other pieces of meta-data that we can check includes the number of features, or the bounding box area the shapefile covers:
>>> len(sf)
663
>>> sf.bbox
[-122.515048, 37.652916, -122.327622, 37.863433]
A shapefile's geometry is the collection of points or shapes made from vertices and implied arcs representing physical locations. All types of shapefiles just store points. The metadata about the points determine how they are handled by software.
You can get the a list of the shapefile's geometry by calling the shapes() method.
>>> shapes = sf.shapes()
The shapes method returns a list of Shape objects describing the geometry of each shape record.
>>> len(shapes)
663
Each shape record contains the following attributes:
>>> for name in dir(shapes[3]):
... if not name.startswith('__'):
... name
'bbox'
'parts'
'points'
'shapeType'
-
shapeType: an integer representing the type of shape as defined by the shapefile specification.
>>> shapes[3].shapeType 5
-
bbox: If the shape type contains multiple points this tuple describes the lower left (x,y) coordinate and upper right corner coordinate creating a complete box around the points. If the shapeType is a Null (shapeType == 0) then an AttributeError is raised.
>>> # Get the bounding box of the 4th shape. >>> # Round coordinates to 3 decimal places >>> bbox = shapes[3].bbox >>> ['%.3f' % coord for coord in bbox] ['-122.486', '37.787', '-122.446', '37.811']
-
parts: Parts simply group collections of points into shapes. If the shape record has multiple parts this attribute contains the index of the first point of each part. If there is only one part then a list containing 0 is returned.
>>> shapes[3].parts [0]
-
points: The points attribute contains a list of tuples containing an (x,y) coordinate for each point in the shape.
>>> len(shapes[3].points) 173 >>> # Get the 8th point of the fourth shape >>> # Truncate coordinates to 3 decimal places >>> shape = shapes[3].points[7] >>> ['%.3f' % coord for coord in shape] ['-122.471', '37.787']
To read a single shape by calling its index use the shape() method. The index is the shape's count from 0. So to read the 8th shape record you would use its index which is 7.
>>> s = sf.shape(7)
>>> # Read the bbox of the 8th shape to verify
>>> # Round coordinates to 3 decimal places
>>> ['%.3f' % coord for coord in s.bbox]
['-122.450', '37.801', '-122.442', '37.808']
A record in a shapefile contains the attributes for each shape in the collection of geometry. Records are stored in the dbf file. The link between geometry and attributes is the foundation of all geographic information systems. This critical link is implied by the order of shapes and corresponding records in the shp geometry file and the dbf attribute file.
The field names of a shapefile are available as soon as you read a shapefile. You can call the "fields" attribute of the shapefile as a Python list. Each field is a Python list with the following information:
- Field name: the name describing the data at this column index.
- Field type: the type of data at this column index. Types can be: Character, Numbers, Longs, Dates, or Memo. The "Memo" type has no meaning within a GIS and is part of the xbase spec instead.
- Field length: the length of the data found at this column index. Older GIS software may truncate this length to 8 or 11 characters for "Character" fields.
- Decimal length: the number of decimal places found in "Number" fields.
To see the fields for the Reader object above (sf) call the "fields" attribute:
>>> fields = sf.fields
>>> assert fields == [("DeletionFlag", "C", 1, 0), ["AREA", "N", 18, 5],
... ["BKG_KEY", "C", 12, 0], ["POP1990", "N", 9, 0], ["POP90_SQMI", "N", 10, 1],
... ["HOUSEHOLDS", "N", 9, 0],
... ["MALES", "N", 9, 0], ["FEMALES", "N", 9, 0], ["WHITE", "N", 9, 0],
... ["BLACK", "N", 8, 0], ["AMERI_ES", "N", 7, 0], ["ASIAN_PI", "N", 8, 0],
... ["OTHER", "N", 8, 0], ["HISPANIC", "N", 8, 0], ["AGE_UNDER5", "N", 8, 0],
... ["AGE_5_17", "N", 8, 0], ["AGE_18_29", "N", 8, 0], ["AGE_30_49", "N", 8, 0],
... ["AGE_50_64", "N", 8, 0], ["AGE_65_UP", "N", 8, 0],
... ["NEVERMARRY", "N", 8, 0], ["MARRIED", "N", 9, 0], ["SEPARATED", "N", 7, 0],
... ["WIDOWED", "N", 8, 0], ["DIVORCED", "N", 8, 0], ["HSEHLD_1_M", "N", 8, 0],
... ["HSEHLD_1_F", "N", 8, 0], ["MARHH_CHD", "N", 8, 0],
... ["MARHH_NO_C", "N", 8, 0], ["MHH_CHILD", "N", 7, 0],
... ["FHH_CHILD", "N", 7, 0], ["HSE_UNITS", "N", 9, 0], ["VACANT", "N", 7, 0],
... ["OWNER_OCC", "N", 8, 0], ["RENTER_OCC", "N", 8, 0],
... ["MEDIAN_VAL", "N", 7, 0], ["MEDIANRENT", "N", 4, 0],
... ["UNITS_1DET", "N", 8, 0], ["UNITS_1ATT", "N", 7, 0], ["UNITS2", "N", 7, 0],
... ["UNITS3_9", "N", 8, 0], ["UNITS10_49", "N", 8, 0],
... ["UNITS50_UP", "N", 8, 0], ["MOBILEHOME", "N", 7, 0]]
You can get a list of the shapefile's records by calling the records() method:
>>> records = sf.records()
>>> len(records)
663
Each record is a list containing an attribute corresponding to each field in the field list.
For example in the 4th record of the blockgroups shapefile the 2nd and 3rd fields are the blockgroup id and the 1990 population count of that San Francisco blockgroup:
>>> records[3][1:3]
['060750601001', 4715]
To read a single record call the record() method with the record's index:
>>> rec = sf.record(3)
>>> rec[1:3]
['060750601001', 4715]
You may want to examine both the geometry and the attributes for a record at the same time. The shapeRecord() and shapeRecords() method let you do just that.
Calling the shapeRecords() method will return the geometry and attributes for all shapes as a list of ShapeRecord objects. Each ShapeRecord instance has a "shape" and "record" attribute. The shape attribute is a ShapeRecord object as discussed in the first section "Reading Geometry". The record attribute is a list of field values as demonstrated in the "Reading Records" section.
>>> shapeRecs = sf.shapeRecords()
Let's read the blockgroup key and the population for the 4th blockgroup:
>>> shapeRecs[3].record[1:3]
['060750601001', 4715]
Now let's read the first two points for that same record:
>>> points = shapeRecs[3].shape.points[0:2]
>>> len(points)
2
The shapeRecord() method reads a single shape/record pair at the specified index. To get the 4th shape record from the blockgroups shapefile use the third index:
>>> shapeRec = sf.shapeRecord(3)
The blockgroup key and population count:
>>> shapeRec.record[1:3]
['060750601001', 4715]
>>> points = shapeRec.shape.points[0:2]
>>> len(points)
2
PyShp tries to be as flexible as possible when writing shapefiles while maintaining some degree of automatic validation to make sure you don't accidentally write an invalid file.
PyShp can write just one of the component files such as the shp or dbf file without writing the others. So in addition to being a complete shapefile library, it can also be used as a basic dbf (xbase) library. Dbf files are a common database format which are often useful as a standalone simple database format. And even shp files occasionally have uses as a standalone format. Some web-based GIS systems use an user-uploaded shp file to specify an area of interest. Many precision agriculture chemical field sprayers also use the shp format as a control file for the sprayer system (usually in combination with custom database file formats).
To create a shapefile you add geometry and/or attributes using methods in the Writer class until you are ready to save the file.
Create an instance of the Writer class to begin creating a shapefile:
>>> w = shapefile.Writer()
The shape type defines the type of geometry contained in the shapefile. All of the shapes must match the shape type setting.
There are three ways to set the shape type:
- Set it when creating the class instance.
- Set it by assigning a value to an existing class instance.
- Set it automatically to the type of the first non-null shape by saving the shapefile.
To manually set the shape type for a Writer object when creating the Writer:
>>> w = shapefile.Writer(shapeType=3)
>>> w.shapeType
3
OR you can set it after the Writer is created:
>>> w.shapeType = 1
>>> w.shapeType
1
Because every shape must have a corresponding record it is critical that the number of records equals the number of shapes to create a valid shapefile. You must take care to add records and shapes in the same order so that the record data lines up with the geometry data. For example:
>>> w = shapefile.Writer(shapeType=shapefile.POINT)
>>> w.field("field1", "C")
>>> w.field("field2", "C")
>>> w.record("row", "one")
>>> w.point(1, 1)
>>> w.record("row", "two")
>>> w.point(2, 2)
To help prevent accidental misalignment pyshp has an "auto balance" feature to make sure when you add either a shape or a record the two sides of the equation line up. This way if you forget to update an entry the shapefile will still be valid and handled correctly by most shapefile software. Autobalancing is NOT turned on by default. To activate it set the attribute autoBalance to 1 or True:
>>> w.autoBalance = 1
>>> w.record("row", "three")
>>> w.record("row", "four")
>>> w.point(4, 4)
>>> w.recNum == w.shpNum
True
You also have the option of manually calling the balance() method at any time to ensure the other side is up to date. When balancing is used null shapes are created on the geometry side or records with a value of "NULL" for each field is created on the attribute side. This gives you flexibility in how you build the shapefile. You can create all of the shapes and then create all of the records or vice versa.
>>> w.autoBalance = 0
>>> w.record("row", "five")
>>> w.record("row", "six")
>>> w.record("row", "seven")
>>> w.point(5, 5)
>>> w.point(6, 6)
>>> w.balance()
>>> w.recNum == w.shpNum
True
If you do not use the autobalance or balance method and forget to manually balance the geometry and attributes the shapefile will be viewed as corrupt by most shapefile software.
Geometry is added using one of several convenience methods. The "null" method is used for null shapes, "point" is used for point shapes, "line" for lines, and "poly" is used for polygons and everything else.
Adding a Point shape
Point shapes are added using the "point" method. A point is specified by an x, y value. An optional z (elevation) and m (measure) value can be set if the shapeType is PointZ or PointM.
>>> w = shapefile.Writer()
>>> w.field('name', 'C')
>>> w.point(122, 37)
>>> w.record('point1')
>>> w.save('shapefiles/test/point')
Adding a Polygon shape
Shapefile polygons must have at least 4 points and the last point must be the same as the first. PyShp automatically enforces closed polygons.
>>> w = shapefile.Writer()
>>> w.field('name', 'C')
>>> w.poly(parts=[[[122,37,4,9], [117,36,3,4]], [[115,32,8,8],
... [118,20,6,4], [113,24]]])
>>> w.record('polygon1')
>>> w.save('shapefiles/test/polygon')
Adding a Line shape
A line must have at least two points. Because of the similarities between polygon and line types it is possible to create a line shape using either the "line" or "poly" method.
>>> w = shapefile.Writer()
>>> w.field('name', 'C')
>>> w.line(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
>>> w.poly(parts=[[[1,3],[5,3]]], shapeType=shapefile.POLYLINE)
>>> w.record('line1')
>>> w.record('line2')
>>> w.save('shapefiles/test/line')
Adding a Null shape
Because Null shape types (shape type 0) have no geometry the "null" method is called without any arguments. This type of shapefile is rarely used but it is valid.
>>> w = shapefile.Writer()
>>> w.field('name', 'C')
>>> w.null()
>>> w.record('nullgeom')
>>> w.save('shapefiles/test/null')
Adding from an existing Shape object
Finally, geometry can be added by passing an existing "Shape" object to the "shape" method. This can be particularly useful for copying from one file to another:
>>> r = shapefile.Reader('shapefiles/test/polygon')
>>> w = shapefile.Writer()
>>> w.fields = r.fields[1:] # skip first deletion field
>>> for shaperec in r.iterShapeRecords():
... w.record(*shaperec.record)
... w.shape(shaperec.shape)
>>> w.save('shapefiles/test/copy')
Adding record attributes involves two steps. Step 1 is to create fields to contain attribute values and step 2 is to populate the fields with values for each shape record.
There are several different field types, all of which support storing None values as NULL.
Text fields are created using the 'C' type, and the third 'size' argument can be customized to the expected length of text values to save space:
>>> w = shapefile.Writer()
>>> w.field('TEXT', 'C')
>>> w.field('SHORT_TEXT', 'C', size=5)
>>> w.field('LONG_TEXT', 'C', size=250)
>>> w.null()
>>> w.record('Hello', 'World', 'World'*50)
>>> w.save('shapefiles/test/dtype')
>>> r = shapefile.Reader('shapefiles/test/dtype')
>>> assert r.record(0) == ['Hello', 'World', 'World'*50]
Date fields are created using the 'D' type, and can be created using either date objects, lists, or a YYYYMMDD formatted string. Field length or decimal have no impact on this type:
>>> from datetime import date
>>> w = shapefile.Writer()
>>> w.field('DATE', 'D')
>>> w.null()
>>> w.null()
>>> w.null()
>>> w.null()
>>> w.record(date(1998,1,30))
>>> w.record([1998,1,30])
>>> w.record('19980130')
>>> w.record(None)
>>> w.save('shapefiles/test/dtype')
>>> r = shapefile.Reader('shapefiles/test/dtype')
>>> assert r.record(0) == [date(1998,1,30)]
>>> assert r.record(1) == [date(1998,1,30)]
>>> assert r.record(2) == [date(1998,1,30)]
>>> assert r.record(3) == [None]
Numeric fields are created using the 'N' type (or the 'F' type, which is exactly the same). By default the fourth decimal argument is set to zero, essentially creating an integer field. To store floats you must set the decimal argument to the precision of your choice. To store very large numbers you must increase the field length size to the total number of digits (including comma and minus).
>>> w = shapefile.Writer()
>>> w.field('INT', 'N')
>>> w.field('LOWPREC', 'N', decimal=2)
>>> w.field('MEDPREC', 'N', decimal=10)
>>> w.field('HIGHPREC', 'N', decimal=30)
>>> w.field('FTYPE', 'F', decimal=10)
>>> w.field('LARGENR', 'N', 101)
>>> nr = 1.3217328
>>> w.null()
>>> w.null()
>>> w.record(INT=int(nr), LOWPREC=nr, MEDPREC=nr, HIGHPREC=-3.2302e-25, FTYPE=nr, LARGENR=int(nr)*10**100)
>>> w.record(None, None, None, None, None, None)
>>> w.save('shapefiles/test/dtype')
>>> r = shapefile.Reader('shapefiles/test/dtype')
>>> assert r.record(0) == [1, 1.32, 1.3217328, -3.2302e-25, 1.3217328, 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000]
>>> assert r.record(1) == [None, None, None, None, None, None]
Finally, we can create boolean fields by setting the type to 'L'. This field can take True or False values, or 1 (True) or 0 (False). None is interpreted as missing.
>>> w = shapefile.Writer()
>>> w.field('BOOLEAN', 'L')
>>> w.null()
>>> w.null()
>>> w.null()
>>> w.null()
>>> w.null()
>>> w.null()
>>> w.record(True)
>>> w.record(1)
>>> w.record(False)
>>> w.record(0)
>>> w.record(None)
>>> w.record("Nonesense")
>>> w.save('shapefiles/test/dtype')
>>> r = shapefile.Reader('shapefiles/test/dtype')
>>> r.record(0)
[True]
>>> r.record(1)
[True]
>>> r.record(2)
[False]
>>> r.record(3)
[False]
>>> r.record(4)
[None]
>>> r.record(5)
[None]
You can also add attributes using keyword arguments where the keys are field names.
>>> w = shapefile.Writer()
>>> w.field('FIRST_FLD','C','40')
>>> w.field('SECOND_FLD','C','40')
>>> w.record('First', 'Line')
>>> w.record(FIRST_FLD='First', SECOND_FLD='Line')
File extensions are optional when reading or writing shapefiles. If you specify them PyShp ignores them anyway. When you save files you can specify a base file name that is used for all three file types. Or you can specify a name for one or more file types. In that case, any file types not assigned will not save and only file types with file names will be saved. If you do not specify any file names (i.e. save()), then a unique file name is generated with the prefix "shapefile_" followed by random characters which is used for all three files. The unique file name is returned as a string.
>>> targetName = w.save()
>>> assert("shapefile_" in targetName)
Just as you can read shapefiles from python file-like objects you can also write them.
>>> try:
... from StringIO import StringIO
... except ImportError:
... from io import BytesIO as StringIO
>>> shp = StringIO()
>>> shx = StringIO()
>>> dbf = StringIO()
>>> w.saveShp(shp)
>>> w.saveShx(shx)
>>> w.saveDbf(dbf)
>>> # Normally you would call the "StringIO.getvalue()" method on these objects.
>>> shp = shx = dbf = None
Despite being a lightweight library, PyShp is designed to be able to read and write shapefiles of any size, allowing you to work with hundreds of thousands or even millions of records and complex geometries.
When first creating the Reader class, the library only reads the header information and leaves the rest of the file contents alone. Once you call the records() and shapes() methods however, it will attempt to read the entire file into memory at once. For very large files this can result in MemoryError. So when working with large files it is recommended to use instead the iterShapes(), iterRecords(), or iterShapeRecords() methods instead. These iterate through the file contents one at a time, enabling you to loop through them while keeping memory usage at a minimum.
>>> for shape in sf.iterShapes():
... # do something here
... pass
>>> for rec in sf.iterRecords():
... # do something here
... pass
>>> for shapeRec in sf.iterShapeRecords():
... # do something here
... pass
The shapefile Writer class uses a similar streaming approach to keep memory usage at a minimum, except you don't have change any of your code. The library takes care of this under-the-hood by creating a set of temporary files and immediately writing each geometry and record to disk the moment they are added using shape() or record(). You still have to call save() as usual in order to specify the final location of the output file and in order for the header information to be calculated and written to the beginning of the file.
This means that as long as you are able to iterate through a source file without having to load everything into memory, such as a large CSV table or a large shapefile, you can process and write any number of items, and even merging many different source files into a single large shapefile. If you need to edit or undo any of your writing you would have to read the file back in, one record at a time, make your changes, and write it back out.
The Python __geo_interface__ convention provides a data interchange interface among geospatial Python libraries. The interface returns data as GeoJSON which gives you nice compatibility with other libraries and tools including Shapely, Fiona, and PostGIS. More information on the __geo_interface__ protocol can be found at: https://gist.github.com/sgillies/2217756. More information on GeoJSON is available at http://geojson.org.
>>> s = sf.shape(0)
>>> s.__geo_interface__["type"]
'MultiPolygon'
Just as the library can expose its objects to other applications through the geo interface, it also supports receiving objects with the geo interface from other applications. To write shapes based on GeoJSON objects, simply send an object with the geo interface or a GeoJSON dictionary to the shape() method instead of a Shape object. Alternatively, you can construct a Shape object from GeoJSON using the "geojson_as_shape()" function.
>>> w = shapefile.Writer()
>>> w.field('name', 'C')
>>> w.shape( {"type":"Point", "coordinates":[1,1]} )
>>> w.record('one')
>>> shape = shapefile.geojson_to_shape( {"type":"Point", "coordinates":[2,2]} )
>>> shape.shapeType
1
>>> shape.points
[[2, 2]]
>>> w.shape(shape)
>>> w.record('two')
>>> w.save('shapefiles/test/geojson')
The testing framework is doctest, which are located in this file README.md. In the same folder as README.md and shapefile.py, from the command line run
$ python shapefile.py
Linux/Mac and similar platforms will need to run $ dos2unix README.md
in order
correct line endings in README.md.