Change shapefile coordinate system using Python

Question

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)

Answer 1

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

Answer 2

Here is even more concise option using geopandas library:

import geopandas

data = geopandas.read_file(path+file_in)
# change CRS to epsg 4326
data = data.to_crs(epsg=4326)
# write shp file
data.to_file(path+file_out)

BTW Kartograph is not developed anymore according to developer git-hub page!

Answer 3

By using fiona and shapely the code is more concise:

from kartograph import Kartograph
from pyproj import Proj, transform
import fiona
from shapely.geometry import mapping, Point

from fiona.crs import from_epsg
crs = from_epsg(4326)

K = Kartograph()
# S-JTSK EPSG:5514
# Kartograph uses EPSG:4326
inProj = Proj(init='epsg:32612', preserve_units=True)
outProj = Proj(init='epsg:4326')

path = '/home/zeito/pyqgis_data/random_points.shp'

points = fiona.open(path)

new_points = []

for point in points:
    x, y = point['geometry']['coordinates']
    x2,y2 = transform(inProj, outProj, x, y)
    new_points.append((x2, y2))

# creation of the resulting shapefile
schema1 = {'geometry': 'Point','properties': {'id': 'int'},}
with fiona.open('/home/zeito/pyqgis_data/reprojected_points.shp', 'w', 'ESRI Shapefile', schema1, crs)  as output:

    for i, point in enumerate(new_points):
        output.write({'geometry':mapping(Point(point[0], point[1])),'properties': {'id':i}})

I tried it out with point shapefile indicated at code as random_points.shp (EPSG:32612). After running it at Python Console, I loaded reprojected shapefile at Map View of QGIS and it was transformed adequately (as it can be observed at following image).

Now, it is available to use with kartograph python module.

Answer 4

I was having the same problem with the code in this link; https://pcjericks.github.io/py-gdalogr-cookbook/projection.html#reproject-a-layer.

Researching I found the solution here; https://github.com/OSGeo/gdal/issues/1546#issuecomment-616344351

then the solution is to add the following line;

inSpatialRef.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)