使用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