How to use a specific transformation for a coordinate projection with Python

Question

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 03 '19 at 12:16
You would like to transform from DHDN (EPSG:4314) to WGS84 (EPSG:4326) -? — Jane, Jul 03 '19 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 03 '19 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 http://proj.maptools.org/gen_parms.html — mkennedy, Jul 03 '19 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 04 '19 at 08:41
The problem is, that the projected coordinates are nearly correct but not exactly. This is where the transformation-method (NTv2) is needed. I tried to unterstand the comment and the link of mkennedy and I think I understood now, that I need to use the cs2cs tool. But at the moment I don't understand how I can put this into python. — AnPy, Jul 04 '19 at 08:58
@mkennedy: I'm afraid I don't understand exactly what you mean. I downloaded the ".gsb"link file you are talking about. But how can I use this in my script? — AnPy, Jul 04 '19 at 09:19
see here for an example (note - not tested personally): https://jvmora.wordpress.com/2013/02/22/ed50-to-etrs89-using-a-spanish-shift-grid-with-pyproj/ — JimT, Jul 04 '19 at 12:24
This looks good at the moment. Tomorrow morning I will check the results and let you know, if the result ist correct and precisly. Thanks a lot so far. — AnPy, Jul 04 '19 at 15:15

score 3 · Answer 1 · answered Jul 18 '19 at 13:23

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

Please mind that FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string). E.g. it must be now output_proj = pyproj.Proj("EPSG:"+output_EPSG) — Taras, Jan 13 '22 at 12:33

How to use a specific transformation for a coordinate projection with Python

1 Answers1