Save points to a shapefile with GDAL

Save point coordinates and their spatial reference to a shapefile. This tutorial was first presented as part of a three-hour session on Working with Geographic Information Systems in Python during the 2009 Python Conference in Chicago, Illinois.

Loading other geometries such as lines and polygons and their associated fields are covered in Save and load geometries and fields from a shapefile with GDAL.

Requirements

• Python
• Geospatial Data Abstraction Library

For the walkthrough, you will also need the following:

• IPython

Example

Download the code and data, unzip and start Python.

wget http://invisibleroads.com/tutorials/_downloads/gdal-shapefile-points.zip
unzip gdal-shapefile-points.zip
cd gdal-shapefile-points
python

Save points and spatial reference.

import point_store
point_store.save('points-shifted.shp', [
    (474595, 4429281),
    (474580, 4429137),
    (474499, 4429249),
], '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

Concepts

• A point is a type of geometry stored as a feature.
• A layer can have many features.
• A datasource can have many layers.
• The driver saves the datasource in a specific format.
• GDAL has drivers for different formats.

Driver
    Datasource
        Layer
            Feature
                Geometry
                    Point

Walkthrough

Import GDAL.

import osgeo.ogr, osgeo.osr

Set spatial reference.

spatialReference = osgeo.osr.SpatialReference()
spatialReference.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

Create shapefile.

driver = osgeo.ogr.GetDriverByName('ESRI Shapefile')
shapeData = driver.CreateDataSource('points-shifted.shp')

Create layer.

layer = shapeData.CreateLayer('layer1', spatialReference, osgeo.ogr.wkbPoint)
layerDefinition = layer.GetLayerDefn()

Create point.

point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
point.SetPoint(0, 474595, 4429281)

Put point as a geometry inside a feature.

featureIndex = 0
feature = osgeo.ogr.Feature(layerDefinition)
feature.SetGeometry(point)
feature.SetFID(featureIndex)

Put feature in a layer.

layer.CreateFeature(feature)

Flush.

shapeData.Destroy()

You can peek at the contents of the shapefile using Quantum GIS.

Code

'Save and load points to a shapefile'
# Import system modules
import osgeo.ogr
import osgeo.osr
import os


# Core

def save(shapePath, geoLocations, proj4):
    'Save points in the given shapePath'
    # Get driver
    driver = osgeo.ogr.GetDriverByName('ESRI Shapefile')
    # Create shapeData
    shapePath = validateShapePath(shapePath)
    if os.path.exists(shapePath): 
        os.remove(shapePath)
    shapeData = driver.CreateDataSource(shapePath)
    # Create spatialReference
    spatialReference = getSpatialReferenceFromProj4(proj4)
    # Create layer
    layerName = os.path.splitext(os.path.split(shapePath)[1])[0]
    layer = shapeData.CreateLayer(layerName, spatialReference, osgeo.ogr.wkbPoint)
    layerDefinition = layer.GetLayerDefn()
    # For each point,
    for pointIndex, geoLocation in enumerate(geoLocations):
        # Create point
        geometry = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
        geometry.SetPoint(0, geoLocation[0], geoLocation[1])
        # Create feature
        feature = osgeo.ogr.Feature(layerDefinition)
        feature.SetGeometry(geometry)
        feature.SetFID(pointIndex)
        # Save feature
        layer.CreateFeature(feature)
        # Cleanup
        geometry.Destroy()
        feature.Destroy()
    # Cleanup
    shapeData.Destroy()
    # Return
    return shapePath

def load(shapePath):
    'Given a shapePath, return a list of points in GIS coordinates'
    # Open shapeData
    shapeData = osgeo.ogr.Open(validateShapePath(shapePath))
    # Validate shapeData
    validateShapeData(shapeData)
    # Get the first layer
    layer = shapeData.GetLayer()
    # Initialize
    points = []
    # For each point,
    for index in xrange(layer.GetFeatureCount()):
        # Get
        feature = layer.GetFeature(index)
        geometry = feature.GetGeometryRef()
        # Make sure that it is a point
        if geometry.GetGeometryType() != osgeo.ogr.wkbPoint: 
            raise ShapeDataError('This module can only load points; use geometry_store.py')
        # Get pointCoordinates
        pointCoordinates = geometry.GetX(), geometry.GetY()
        # Append
        points.append(pointCoordinates)
        # Cleanup
        feature.Destroy()
    # Get spatial reference as proj4
    proj4 = layer.GetSpatialRef().ExportToProj4()
    # Cleanup
    shapeData.Destroy()
    # Return
    return points, proj4

def merge(sourcePaths, targetPath):
    'Merge a list of shapefiles into a single shapefile'
    # Load
    items = [load(validateShapePath(x)) for x in sourcePaths]
    pointLists = [x[0] for x in items]
    points = reduce(lambda x,y: x+y, pointLists)
    spatialReferences= [x[1] for x in items]
    # Make sure that all the spatial references are the same
    if len(set(spatialReferences)) != 1: 
        raise ShapeDataError('The shapefiles must have the same spatial reference')
    spatialReference = spatialReferences[0]
    # Save
    save(validateShapePath(targetPath), points, spatialReference)

def getSpatialReferenceFromProj4(proj4):
    'Return GDAL spatial reference object from proj4 string'
    spatialReference = osgeo.osr.SpatialReference()
    spatialReference.ImportFromProj4(proj4)
    return spatialReference


# Validate

def validateShapePath(shapePath):
    'Validate shapefile extension'
    return os.path.splitext(str(shapePath))[0] + '.shp'

def validateShapeData(shapeData):
    'Make sure we can access the shapefile'
    # Make sure the shapefile exists
    if not shapeData:
        raise ShapeDataError('The shapefile is invalid')
    # Make sure there is exactly one layer
    if shapeData.GetLayerCount() != 1:
        raise ShapeDataError('The shapefile must have exactly one layer')


# Error

class ShapeDataError(Exception):
    pass