投影变换

标签：

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os

try:
import ogr
import osr

except ImportError:
from osgeo import ogr, osr

shp_driver = ogr.GetDriverByName(‘ESRI Shapefile‘)

# input SpatialReference
input_srs = osr.SpatialReference()
input_srs.ImportFromEPSG(4326)

# output SpatialReference
output_srs = osr.SpatialReference()
output_srs.ImportFromEPSG(3857)

# create the CoordinateTransformation
coord_trans = osr.CoordinateTransformation(input_srs, output_srs)

# get the input layer
input_shp = shp_driver.Open(r‘../geodata/UTM_Zone_Boundaries.shp‘)
in_shp_layer = input_shp.GetLayer()

# create the output layer
output_shp_file = r‘../geodata/UTM_Zone_Boundaries_3857.shp‘
# check if output file exists if yes delete it
if os.path.exists(output_shp_file):
shp_driver.DeleteDataSource(output_shp_file)

# create a new Shapefile object
output_shp_dataset = shp_driver.CreateDataSource(output_shp_file)

# create a new layer in output Shapefile and define its geometry type
output_shp_layer = output_shp_dataset.CreateLayer("UTM_Zone_Boundaries_3857", geom_type=ogr.wkbMultiPolygon)

# add fields to the new output Shapefile
# get list of attribute fields
# create new fields for output
in_layer_def = in_shp_layer.GetLayerDefn()
for i in range(0, in_layer_def.GetFieldCount()):
field_def = in_layer_def.GetFieldDefn(i)
output_shp_layer.CreateField(field_def)

# get the output layer‘s feature definition
output_layer_def = output_shp_layer.GetLayerDefn()

# loop through the input features
in_feature = in_shp_layer.GetNextFeature()
while in_feature:
# get the input geometry
geom = in_feature.GetGeometryRef()
# reproject the geometry
geom.Transform(coord_trans)
# create a new feature
output_feature = ogr.Feature(output_layer_def)
# set the geometry and attribute
output_feature.SetGeometry(geom)
for i in range(0, output_layer_def.GetFieldCount()):
output_feature.SetField(output_layer_def.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))
# add the feature to the shapefile
output_shp_layer.CreateFeature(output_feature)
# destroy the features and get the next input feature
output_feature.Destroy()
in_feature.Destroy()
in_feature = in_shp_layer.GetNextFeature()

# close the shapefiles
input_shp.Destroy()
output_shp_dataset.Destroy()

spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(3857)

spatialRef.MorphToESRI()
prj_file = open(‘UTM_Zone_Boundaries.prj‘, ‘w‘)
prj_file.write(spatialRef.ExportToWkt())
prj_file.close()

投影变换

标签：

原文地址：http://www.cnblogs.com/gispathfinder/p/5722103.html