#!/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()