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)
  • Please Edit your code to be legible by using the Code formatter ({}) to indent all lines an additional four spaces. The backticks are only used inline with normal text. – Vince Dec 16 '17 at 15:06
  • @vesimies.h You can use fiona and shapely. Code is more concise. – xunilk Dec 17 '17 at 18:03

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()

# output SpatialReference
outSpatialRef = osr.SpatialReference()

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

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

# create the output layer
if os.path.exists(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)

# 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
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    # dereference the features and get the next input feature
    outFeature = None
    inFeature = inLayer.GetNextFeature()

# Save and close the shapefiles
inDataSet = None
outDataSet = None
  • This doesn't work in gdal version 3.0.2. I just found that out the hard way. It does seem to work in gdal version 2.3.3. – C Johnson Feb 20 at 17:12

Here is even more concise option using geopandas library:

import geopandas

data = geopandas.read_file(path+file_in)
# change CRS to epsg 4326
## here it is theoretically possible to use fiona.crs function from_epsg(4326), 
## however it doesnt work properly on my Windows PC, so below is proj4 string
data.crs = '+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs' 
# write shp file

It is projected properly, however for some reason it is not recognised in QGIS as epsg 4326 or any official projection, if someone knows the reason I would appreciate your comments:)

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

  • I like this code except that the proj4 string I have for 28992 is not working as I expect. The generated shapefile in arcgis properties shows that it is stereographic projection. However, it doesn't confirm to all standards of 28992 :( Any suggestions?spatialreference.org/ref/epsg/amersfoort-rd-new/proj4 – April Oct 14 '19 at 7:30

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).

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.