使用python对shapefile重投影
脚本如下：
#!/usr/bin/env python
from osgeo import ogr, osr
from osgeo import gdal
import os
def reproject(inputfile,outputfile,layername):
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
gdal.SetConfigOption("SHAPE_ENCODING","")
driver = ogr.GetDriverByName('ESRI Shapefile')
# input SpatialReference
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(4326)
# output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(3857)
# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
# get the input layer
inDataSet = driver.Open(inputfile)
inLayer = inDataSet.GetLayer()
# create the output layer
outputShapefile = outputfile
outDataSet = driver.CreateDataSource(outputShapefile)
print inLayer.GetGeomType()
outLayer = outDataSet.CreateLayer(layername,geom_type=inLayer.GetGeomType() )
# add fields
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outLayer.CreateField(fieldDefn)
# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()
# loop through the input features
inFeature = inLayer.GetNextFeature()
while inFeature:
# get the input geometry
geom = inFeature.GetGeometryRef()
# reproject the geometry
geom.Transform(coordTrans)
# create a new feature
outFeature = ogr.Feature(outLayerDefn)
# set the geometry and attribute
outFeature.SetGeometry(geom)
for i in range(0, outLayerDefn.GetFieldCount()):
outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
# add the feature to the shapefile
outLayer.CreateFeature(outFeature)
# destroy the features and get the next input feature
outFeature.Destroy()
inFeature.Destroy()
inFeature = inLayer.GetNextFeature()
# close the shapefiles
inDataSet.Destroy()
outDataSet.Destroy()
查看原文：http://www.giser.net/?p=1335
转载自：https://blog.csdn.net/barry114/article/details/46648301
