'Save and load points to a shapefile'# Import system modulesimportosgeo.ogrimportosgeo.osrimportos# Coredefsave(shapePath,geoLocations,proj4):'Save points in the given shapePath'# Get driverdriver=osgeo.ogr.GetDriverByName('ESRI Shapefile')# Create shapeDatashapePath=validateShapePath(shapePath)ifos.path.exists(shapePath):os.remove(shapePath)shapeData=driver.CreateDataSource(shapePath)# Create spatialReferencespatialReference=getSpatialReferenceFromProj4(proj4)# Create layerlayerName=os.path.splitext(os.path.split(shapePath)[1])[0]layer=shapeData.CreateLayer(layerName,spatialReference,osgeo.ogr.wkbPoint)layerDefinition=layer.GetLayerDefn()# For each point,forpointIndex,geoLocationinenumerate(geoLocations):# Create pointgeometry=osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)geometry.SetPoint(0,geoLocation[0],geoLocation[1])# Create featurefeature=osgeo.ogr.Feature(layerDefinition)feature.SetGeometry(geometry)feature.SetFID(pointIndex)# Save featurelayer.CreateFeature(feature)# Cleanupgeometry.Destroy()feature.Destroy()# CleanupshapeData.Destroy()# ReturnreturnshapePathdefload(shapePath):'Given a shapePath, return a list of points in GIS coordinates'# Open shapeDatashapeData=osgeo.ogr.Open(validateShapePath(shapePath))# Validate shapeDatavalidateShapeData(shapeData)# Get the first layerlayer=shapeData.GetLayer()# Initializepoints=[]# For each point,forindexinxrange(layer.GetFeatureCount()):# Getfeature=layer.GetFeature(index)geometry=feature.GetGeometryRef()# Make sure that it is a pointifgeometry.GetGeometryType()!=osgeo.ogr.wkbPoint:raiseShapeDataError('This module can only load points; use geometry_store.py')# Get pointCoordinatespointCoordinates=geometry.GetX(),geometry.GetY()# Appendpoints.append(pointCoordinates)# Cleanupfeature.Destroy()# Get spatial reference as proj4proj4=layer.GetSpatialRef().ExportToProj4()# CleanupshapeData.Destroy()# Returnreturnpoints,proj4defmerge(sourcePaths,targetPath):'Merge a list of shapefiles into a single shapefile'# Loaditems=[load(validateShapePath(x))forxinsourcePaths]pointLists=[x[0]forxinitems]points=reduce(lambdax,y:x+y,pointLists)spatialReferences=[x[1]forxinitems]# Make sure that all the spatial references are the sameiflen(set(spatialReferences))!=1:raiseShapeDataError('The shapefiles must have the same spatial reference')spatialReference=spatialReferences[0]# Savesave(validateShapePath(targetPath),points,spatialReference)defgetSpatialReferenceFromProj4(proj4):'Return GDAL spatial reference object from proj4 string'spatialReference=osgeo.osr.SpatialReference()spatialReference.ImportFromProj4(proj4)returnspatialReference# ValidatedefvalidateShapePath(shapePath):'Validate shapefile extension'returnos.path.splitext(str(shapePath))[0]+'.shp'defvalidateShapeData(shapeData):'Make sure we can access the shapefile'# Make sure the shapefile existsifnotshapeData:raiseShapeDataError('The shapefile is invalid')# Make sure there is exactly one layerifshapeData.GetLayerCount()!=1:raiseShapeDataError('The shapefile must have exactly one layer')# ErrorclassShapeDataError(Exception):pass