[GIS] Change shapefile coordinate system using Python

coordinate systemkartographpyprojpythonshapefile

I am struggling with Python and Kartograph.py. I'd like to change coordinate system of the whole shapefile from EPSG:5514 to EPSG:4326. Found some code in here working for points, but I don't know how to cope with shapefile. Code for point goes like this:

from kartograph import Kartograph
from pyproj import Proj, transform

K = Kartograph()
# S-JTSK EPSG:5514
# Kartograph uses EPSG:4326

inProj = Proj(init='epsg:5514', preserve_units=True)
outProj = Proj(init='epsg:4326')

x1,y1 = -599411.949672, -599411.949672
x2,y2 = transform(inProj,outProj,x1,y1)

Do I need to do anything like this for shp?

from kartograph import Kartograph
from pyproj import Proj, transform

K = Kartograph()
# S-JTSK EPSG:5514
# Kartograph uses EPSG:4326

inProj = Proj(init='epsg:5514', preserve_units=True)
outProj = Proj(init='epsg:4326')

for point in shapefile:
    SX = #coordinateX in shp attribute table
    SY = #coordinateY
    x1,y1 = -599411.949672, -599411.949672 #do something to transform every point
    x2,y2 = transform(inProj,outProj,x1,y1)

Best Answer

Here's some handy code from the Python GDAL/OGR cookbook that will reproject a shapefile.

from osgeo import ogr, osr
import os

in_epsg = 5514
out_epsg = 4326
in_shp = '/path/to/input.shp'
out_shp = '/path/to/reprojected.shp'

driver = ogr.GetDriverByName('ESRI Shapefile')

# input SpatialReference
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(in_epsg)

# output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(out_epsg)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# get the input layer
inDataSet = driver.Open(in_shp)
inLayer = inDataSet.GetLayer()

# create the output layer
if os.path.exists(out_shp):
    driver.DeleteDataSource(out_shp)
outDataSet = driver.CreateDataSource(out_shp)
outLayer = outDataSet.CreateLayer("reproject", 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)
    # dereference the features and get the next input feature
    outFeature = None
    inFeature = inLayer.GetNextFeature()

# Save and close the shapefiles
inDataSet = None
outDataSet = None
Related Question