3

I need to project EPSG:3068 (Soldner Berlin) to lat/long.

It works fine with https://epsg.io/3068/map, but I'm looking for a way to do it with a Python script.

I tried both pyproj and the epsg.io API:

from pyproj import Proj
import json
import urllib2

# projection with pyproj
p1 = Proj(init='epsg:3068')
lon, lat = p1(26912,16806,inverse=True)

# projection with the espg.io API
data = json.load(urllib2.urlopen('http://epsg.io/trans?x=26912&y=16608&z=0&s_srs=3068&t_srs=4326'))

But both are a bit off the actual location (in Berlin): enter image description here

The epsg.io marker in the middle is correct. The upper right one is the location that I get from pyproj. The lower left marker is from the epsg.io API.

How can I avoid the discrepancies and get the correct result in Python?

P.S.: Code for creating the map:

import folium
map_1 = folium.Map(location=[lat, lon], zoom_start=100)
map_1.simple_marker([lat, lon], popup='pyproj')
map_1.simple_marker([52.47825965, 13.4327974], popup='epsg.io')
map_1.simple_marker([float(data['y']),float(data['x'])], popup='epsg.io API')
map_1.create_map(path='mthood.html')

EDIT

This solution based on the answer below works:

p1 = Proj(init='epsg:3068')
p2 = Proj(init='epsg:4326')
lon, lat = transform(p1,p2, 26912,16806)
  • 1
    Looks like a typo in the URL (or pyproj): y=16608 or lon, lat = p1(26912,16806,inverse=True) – BradHards Oct 14 '15 at 10:46
  • Oops, indeed! Thanks @BradHards. So the API location is correct now. – tobip Oct 14 '15 at 19:27
2

In your pyproj code, you have not requested the result to be in EPSG:4326 latlong. So Pyproj just transforms Soldner X/Y to latlong on the same ellipsoid (DHDN Bessel, EPSG:4314). The result will be offset against the WGS84 ellipsoid unless you transform it to EPSG:4326.

You can use the proper proj.4 definition of DHDN Soldner Berlin:

+proj=cass +lat_0=52.41864827777778 +lon_0=13.62720366666667 +x_0=40000 +y_0=10000 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs
  • Works! I added the final code above. – tobip Oct 14 '15 at 19:48
1

The difference can be cause by the different SRS definitions for the different environments. You should check the EPSG definition of the projection, for example in my proj4 install in the epsg file it looks like this:

<3068> +proj=cass +lat_0=52.41864827777778 +lon_0=13.62720366666667 +x_0=40000 +y_0=10000 +ellps=bessel +datum=potsdam +units=m +no_defs <>

In QGIS I have the following:

+proj=cass +lat_0=52.41864827777778 +lon_0=13.62720366666667 +x_0=40000 +y_0=10000 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs

May be the different programming environments maintain different epsg definition files.

You should look after the +towgs84 parameter (7 parameter 3D transformation) what can be used to make the transformation more precise.

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.