[GIS] How to use a specific transformation for a coordinate projection with Python

coordinate systemprojpython

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

Best Answer

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

Related Question