6

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

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
3
  • 1
    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 '20 at 17:12
  • How does it not work? Oct 17 '20 at 8:35
  • I just made a tiny modification: outLayer = outDataSet.CreateLayer("reproject", geom_type=ogr.wkbMultiPolygon) for outLayer = outDataSet.CreateLayer( "reproject", outSpatialRef, geom_type=ogr.wkbMultiPolygon). It was to specify the CRS in output file.
    – Rafa
    Aug 25 '21 at 21:03
6

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!

2
  • 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
  • you project the data in wgs84, but then you use a proj string to define the projection as a UTM projection. Therefor your coordinate value are coorect, but the crs is incorrect.
    – radouxju
    Feb 25 '21 at 13:34
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).

enter image description here

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

-2

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)

1
  • 1
    If you have a NEW question, please ask it by clicking the Ask Question button. If you have sufficient reputation, you may upvote the question. Alternatively, "star" it as a favorite and you will be notified of any new answers.
    – Ian Turton
    May 16 '20 at 10:08

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.