2

I try to use Python to project single coordinates from Gauß-Krueger-coordinate system (EPSG:31467) to ETRS89-UTM (EPSG:25832). I'm using the way Antonio Falciano showed here.

So my script looks like this:

from pyproj import Proj, transform

x_in = 3468511.6749226563
y_in = 5428925.356241324

input_proj = Proj(init='epsg:31467')
output_proj = Proj(init='epsg:25832')

x_out, y_out = transform(input_proj, output_proj, x_in, y_in)
print x_out, y_out

returns 468448.880312 5427193.93989

This works for me, but now I want to add a specific geographic (datum) transfomation (DHDN_To_WGS_1984_4_NTv2). How can I manage this. I would prefere to do this without using a GIS but I cannot find a way so far.

After the tip of mkenedy I found out, that I have to use cs2cs and a .gsb-file. I think the following line has to put in my code somehow, but I really have no idea how. If I try, it'll cause a SyntaxError: invalid syntax:

cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel + units=m + nadgrids=BETA2007.gsb +to +proj=utm +ellps=GRS80 +zone=32 +nadgrids=@null x_in x_out 

Although i have been reading internet pages about cs2cs for hours now, I don't understand how to use cs2cs correctly. Can someone help me with this?

Now I tried to do an reprojection first like it is done on this site. It looks like this:

from pyproj import Proj, transform
import pyproj

# GK3 = EPSG: 31467
input_EPSG = '31467'
# UTM32 = EPSG: 25832
output_EPSG = '25832'
gsb_filename = 'BWTA2017.gsb'
x_in = 3468511.6749226563
y_in = 5428925.356241324

gk = pyproj.Proj("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +nadgrids=" + gsb_filename + "+units=m + no_defs")
x, y = gk(x_in, y_in)
print x, y

returns: 1e+30 1e+30
RuntimeError: latitude or longitude exceeded limits
  • Explain please what this specific transfomation means? I found this is a Geographic(datum) transformation, do you have proj4 definition of this or anything else? – Jane Jul 3 at 12:16
  • You would like to transform from DHDN (EPSG:4314) to WGS84 (EPSG:4326) -? – Jane Jul 3 at 12:22
  • No, sorry, I apologize for the confusion. The input system is the DHDN (EPSG: 31467) and it has to be in ETRS89-UTM (EPSG: 25832) – AnPy Jul 3 at 14:13
  • I don't know for sure, but I think you would have to change your Proj definitions to use the proj4 strings and include +nadgrids=dhdnfilename.gsb parameter which means you have to have that grid file downloaded and available. Partial info at proj.maptools.org/gen_parms.html – mkennedy Jul 3 at 16:57
  • @AnPy, EPSG: 31467 and EPSG: 25832 are in your code, and you mentioned "This works for you", So what specific transformation you would like to add? – Jane Jul 4 at 8:41
1

!!!Solved!!! The way that worked in the end:

import pyproj

x_in = 3468072.147662042
y_in = 5430007.049974753

input_EPSG = "31467"
output_EPSG = "25832"

def proj_coord(x_in, y_in):
    """
    projection process - In this method the projection of a single pair of coordinates takes place.
    :param x_in: X - coordinate-value
    :param y_in: Y - coordinate-value
    :return: x_out, y_out (projected and transformed coordinate)
    """
    input_proj = pyproj.Proj(init="EPSG:"+input_EPSG, nadgrids='..\\' + gsb_filename)
    output_proj = pyproj.Proj(init="EPSG:"+output_EPSG)
    print "Coordinates are currently in EPSG ", input_EPSG

    x_out, y_out = pyproj.transform(input_proj, output_proj, x_in, y_in)

    print "Now the coordinates are projected to EPSG", output_EPSG
    print "\t The new coordinate - X: ", x_out
    print "\t The new coordinate - Y: ", y_out
    return x_out, y_out

returns: 
Coordinates are currently in EPSG  31467
Now the coordinates are projected to EPSG 25832
    The new coordinate - X:  468009.548038
    The new coordinate - Y:  5428274.49364

Thanks to all of you. Your questions and tips helped bring to find the right ideas and this procedural method.

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.