from osgeo import ogr, osr import os driver = ogr.GetDriverByName('ESRI Shapefile') # input SpatialReference inSpatialRef = osr.SpatialReference() inSpatialRef.ImportFromEPSG(2927) # output SpatialReference outSpatialRef = osr.SpatialReference() outSpatialRef.ImportFromEPSG(4326) # create the CoordinateTransformation coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef) # get the input layer inDataSet = driver.Open(r'c:\data\spatial\basemap.shp') inLayer = inDataSet.GetLayer() # create the output layer outputShapefile = r'c:\data\spatial\basemap_4326.shp' if os.path.exists(outputShapefile): driver.DeleteDataSource(outputShapefile) outDataSet = driver.CreateDataSource(outputShapefile) outLayer = outDataSet.CreateLayer("basemap_4326", geom_type=ogr.wkbMultiPolygon) # 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()
投影的一层
关注公众号
获取免费资源
Copyright © Since 2014.
开源地理空间基金会中文分会
吉ICP备05002032号
Powered by TorCMS