• Willkommen im Geoclub - dem größten deutschsprachigen Geocaching-Forum. Registriere dich kostenlos, um alle Inhalte zu sehen und neue Beiträge zu erstellen.

Gauß-Krüger in SQL Datenbank konvertieren

Higgins

Geonewbie
Hallo zusammen,

ich stehe vor der Aufgabe, Gauß-Krüger Koordinaten aus einer Datenbanktabelle für die Anzeige mit dem Bing Kartendienst kompatibel zu machen. Dazu muss ich sie konvertieren, am besten wäre wohl Längen- und Breitengrad.

Ich habe dazu im Internet was gefunden, was grundsätzlich auch zu funktionieren scheint. Allerdings ist vermutlich der Bezugsmeridian falsch oder so, ich bekomme nicht die richtigen Ergebnisse.

Hat hier vielleicht jemand eine Idee, welche Werte ich anpassen muss, damit das Ganze für Hamburg funktioniert?

ich vermute stark, dass ich in diesem Bereich was anpassen muss:
, @K0 = 0.9996012717 -- grid scale factor on central meridean
, @OriginLat = 49.0
, @OriginLong = -2.0
, @OriginX = 400000 -- 400 kM
, @OriginY = -100000 -- 100 kM
, @a = 6377563.396 -- Airy Spheroid
, @b = 6356256.910

Und natürlich hier:
, @height = 24 --London's mean height above sea level is 24 metres. Adjust for other locations.

Wäre toll, wenn hier jemand weiterhelfen kann.

Hier der gesamte Code:
****** Object: UserDefinedFunction [dbo].[NEtoLL] Script Date: 09/06/2012 17:06:39 ******/
SET ANSI_NULLS ON
GO

SET QUOTED_IDENTIFIER ON
GO

CREATE FUNCTION [dbo].[NEtoLL] (@East INT, @North INT, @LatOrLng VARCHAR(3)) RETURNS FLOAT AS
BEGIN

--Author: Sandy Motteram
--Date: 06 September 2012

--UDF adapted from javascript at http://www.bdcc.co.uk/LatLngToOSGB.js
--found on page http://mapki.com/wiki/Tools:Snippets

--Instructions:
--Latitude and Longitude are calculated based on BOTH the easting and northing values from the OSGB36
--This UDF takes both easting and northing values in OSGB36 projection and you must specify if a latitude or longitude co-ordinate should be returned.
--IT first converts E/N values to lat and long in OSGB36 projection, then converts those values to lat/lng in WGS84 projection

--Sample values below
--DECLARE @East INT, @North INT, @LatOrLng VARCHAR(3)
--SELECT @East = 529000, @North = 183650 --that combo should be the corner of Camden High St and Delancey St


DECLARE @Pi FLOAT
, @K0 FLOAT
, @OriginLat FLOAT
, @OriginLong FLOAT
, @OriginX FLOAT
, @OriginY FLOAT
, @a FLOAT
, @b FLOAT
, @e2 FLOAT
, @ex FLOAT
, @n1 FLOAT
, @n2 FLOAT
, @n3 FLOAT
, @OriginNorthings FLOAT
, @lat FLOAT
, @lon FLOAT
, @Northing FLOAT
, @Easting FLOAT

SELECT @Pi = 3.14159265358979323846
, @K0 = 0.9996012717 -- grid scale factor on central meridean
, @OriginLat = 49.0
, @OriginLong = -2.0
, @OriginX = 400000 -- 400 kM
, @OriginY = -100000 -- 100 kM
, @a = 6377563.396 -- Airy Spheroid
, @b = 6356256.910
/* , @e2
, @ex
, @n1
, @n2
, @n3
, @OriginNorthings*/

-- compute interim values
SELECT @a = @a * @K0
, @b = @b * @K0

SET @n1 = (@a - @b) / (@a + @b)
SET @n2 = @n1 * @n1
SET @n3 = @n2 * @n1

SET @lat = @OriginLat * @Pi / 180.0 -- to radians

SELECT @e2 = (@a * @a - @b * @b) / (@a * @a) -- first eccentricity
, @ex = (@a * @a - @b * @b) / (@b * @b) -- second eccentricity

SET @OriginNorthings = @b * @lat + @b * (@n1 * (1.0 + 5.0 * @n1 * (1.0 + @n1) / 4.0) * @lat
- 3.0 * @n1 * (1.0 + @n1 * (1.0 + 7.0 * @n1 / 8.0)) * SIN(@lat) * COS(@lat)
+ (15.0 * @n1 * (@n1 + @n2) / 8.0) * SIN(2.0 * @lat) * COS(2.0 * @lat)
- (35.0 * @n3 / 24.0) * SIN(3.0 * @lat) * COS(3.0 * @lat))

SELECT @northing = @north - @OriginY
, @easting = @east - @OriginX

DECLARE @nu FLOAT
, @phid FLOAT
, @phid2 FLOAT
, @t2 FLOAT
, @t FLOAT
, @q2 FLOAT
, @c FLOAT
, @s FLOAT
, @nphid FLOAT
, @dnphid FLOAT
, @nu2 FLOAT
, @nudivrho FLOAT
, @invnurho FLOAT
, @rho FLOAT
, @eta2 FLOAT

/* Evaluate M term: latitude of the northing on the centre meridian */

SET @northing = @northing + @OriginNorthings

SET @phid = @northing / (@b*(1.0 + @n1 + 5.0 * (@n2 + @n3) / 4.0)) - 1.0
SET @phid2 = @phid + 1.0

WHILE (ABS(@phid2 - @phid) > 0.000001)
BEGIN
SET @phid = @phid2;
SET @nphid = @b * @phid + @b * (@n1 * (1.0 + 5.0 * @n1 * (1.0 + @n1) / 4.0) * @phid
- 3.0 * @n1 * (1.0 + @n1 * (1.0 + 7.0 * @n1 / 8.0)) * SIN(@phid) * COS(@phid)
+ (15.0 * @n1 * (@n1 + @n2) / 8.0) * SIN(2.0 * @phid) * COS(2.0 * @phid)
- (35.0 * @n3 / 24.0) * SIN(3.0 * @phid) * COS(3.0 * @phid))

SET @dnphid = @b * ((1.0 + @n1 + 5.0 * (@n2 + @n3) / 4.0) - 3.0 * (@n1 + @n2 + 7.0 * @n3 / 8.0) * COS(2.0 * @phid)
+ (15.0 * (@n2 + @n3) / 4.0) * COS(4 * @phid) - (35.0 * @n3 / 8.0) * COS(6.0 * @phid))

SET @phid2 = @phid - (@nphid - @northing) / @dnphid
END

SELECT @c = COS(@phid)
, @s = SIN(@phid)
, @t = TAN(@phid)
SELECT @t2 = @t * @t
, @q2 = @easting * @easting

SET @nu2 = (@a * @a) / (1.0 - @e2 * @s * @s)
SET @nu = SQRT(@nu2)

SET @nudivrho = @a * @a * @c * @c / (@b * @b) - @c * @c + 1.0
SET @eta2 = @nudivrho - 1
SET @rho = @nu / @nudivrho;

SET @invnurho = ((1.0 - @e2 * @s * @s) * (1.0 - @e2 * @s * @s)) / (@a * @a * (1.0 - @e2))

SET @lat = @phid - @t * @q2 * @invnurho / 2.0 + (@q2 * @q2 * (@t / (24 * @rho * @nu2 * @nu) * (5 + (3 * @t2) + @eta2 - (9 * @t2 * @eta2))))
SET @lon = (@easting / (@c * @nu))
- (@easting * @q2 * ((@nudivrho + 2.0 * @t2) / (6.0 * @nu2)) / (@c * @nu))
+ (@q2 * @q2 * @easting * (5 + (28 * @t2) + (24 * @t2 * @t2)) / (120 * @nu2 * @nu2 * @nu * @c))


SELECT @lat = @lat * 180.0 / @Pi
, @lon = @lon * 180.0 / @Pi + @OriginLong


--Now convert the lat and long from OSGB36 to WGS84

DECLARE @OGlat FLOAT
, @OGlon FLOAT
, @height FLOAT

SELECT @OGlat = @lat
, @OGlon = @lon
, @height = 24 --London's mean height above sea level is 24 metres. Adjust for other locations.

DECLARE @deg2rad FLOAT
, @rad2deg FLOAT
, @radOGlat FLOAT
, @radOGlon FLOAT

SELECT @deg2rad = @Pi / 180
, @rad2deg = 180 / @Pi

--first off convert to radians
SELECT @radOGlat = @OGlat * @deg2rad
, @radOGlon = @OGlon * @deg2rad
--these are the values for WGS84(GRS80) to OSGB36(Airy)

DECLARE @a2 FLOAT
, @h FLOAT
, @xp FLOAT
, @yp FLOAT
, @zp FLOAT
, @xr FLOAT
, @yr FLOAT
, @zr FLOAT
, @sf FLOAT
, @e FLOAT
, @v FLOAT
, @x FLOAT
, @y FLOAT
, @z FLOAT
, @xrot FLOAT
, @yrot FLOAT
, @zrot FLOAT
, @hx FLOAT
, @hy FLOAT
, @hz FLOAT
, @newLon FLOAT
, @newLat FLOAT
, @p FLOAT
, @errvalue FLOAT
, @lat0 FLOAT

SELECT @a2 = 6378137 -- WGS84_AXIS
, @e2 = 0.00669438037928458 -- WGS84_ECCENTRIC
, @h = @height -- height above datum (from $GPGGA sentence)
, @a = 6377563.396 -- OSGB_AXIS
, @e = 0.0066705397616 -- OSGB_ECCENTRIC
, @xp = 446.448
, @yp = -125.157
, @zp = 542.06
, @xr = 0.1502
, @yr = 0.247
, @zr = 0.8421
, @s = -20.4894

-- convert to cartesian; lat, lon are in radians
SET @sf = @s * 0.000001
SET @v = @a / (sqrt(1 - (@e * (SIN(@radOGlat) * SIN(@radOGlat)))))
SET @x = (@v + @h) * COS(@radOGlat) * COS(@radOGlon)
SET @y = (@v + @h) * COS(@radOGlat) * SIN(@radOGlon)
SET @z = ((1 - @e) * @v + @h) * SIN(@radOGlat)

-- transform cartesian
SET @xrot = (@xr / 3600) * @deg2rad
SET @yrot = (@yr / 3600) * @deg2rad
SET @zrot = (@zr / 3600) * @deg2rad
SET @hx = @x + (@x * @sf) - (@y * @zrot) + (@z * @yrot) + @xp
SET @hy = (@x * @zrot) + @y + (@y * @sf) - (@z * @xrot) + @yp
SET @hz = (-1 * @x * @yrot) + (@y * @xrot) + @z + (@z * @sf) + @zp

-- Convert back to lat, lon
SET @newLon = ATAN(@hy / @hx)
SET @p = SQRT((@hx * @hx) + (@hy * @hy))
SET @newLat = ATAN(@hz / (@p * (1 - @e2)))
SET @v = @a2 / (SQRT(1 - @e2 * (SIN(@newLat) * SIN(@newLat))))
SET @errvalue = 1.0;
SET @lat0 = 0
WHILE (@errvalue > 0.001)
BEGIN
SET @lat0 = ATAN((@hz + @e2 * @v * SIN(@newLat)) / @p)
SET @errvalue = ABS(@lat0 - @newLat)
SET @newLat = @lat0
END

--convert back to degrees
SET @newLat = @newLat * @rad2deg
SET @newLon = @newLon * @rad2deg

DECLARE @ReturnMe FLOAT
SET @ReturnMe = 0

IF @LatOrLng = 'Lat'
SET @ReturnMe = @newLat
IF @LatOrLng = 'Lng'
SET @ReturnMe = @newLon

RETURN @ReturnMe
END
GO
 

steffen0815

Geocacher
Hallo,
zum SQL-Server selbst kann ich nix sagen, aber mit einem entsprechenden "räumlichen Aufsatz" auf die Datenbank sollten entsprechende Funktionen auch dort schon vorhanden sein.
Ich denke also nicht, dass man solche Standardfunktionen selbst definieren muss.

Beispiel:
GK Zone 3 (DHDN) : EPSG:31467
Geogr. Lat/lon EPSG:4326
In PostgreSQL/PostGIS wäre es z.B. so:
Code:
select ST_AsText(ST_Transform(st_setsrid(st_makepoint(3570000,5950000),31467),4326))

Die Transformation muss/soll innerhalb der Datenbank erfolgen?


Gruß Steffen
 
OP
H

Higgins

Geonewbie
Ach so, ich dachte UTM wäre das gleiche wie GK und ich muss nur Werte anpassen, um in der richtigen Zone zu landen oder so.

Die Transformation muss in der Datenbank erfolgen, ja. Ich habe einen ETL-Prozess, der rund 400.000 Datensätze aus einer MS_SQL Datenbank holt und mit gewissen Umwandlungen und Verknüpfungen in eine andere MS_SQL Datenbank schreibt. In der Quelle habe ich GK Koordinaten, im Ziel brauche ich ein Format, das ich in Power BI an Bing Maps anbinden kann.

Was ist denn mit "räumlichen Aufsatz" gemeint? Eine zusätzliche Installation oder Import einer Bibliothek oder sowas? Meine Datenbank kennt jedenfalls keine integrierte Funktion namens st_makepoint (und auch die anderen st_funktionsname nicht).
 

steffen0815

Geocacher
Hallo,
ich habe noch mal etwas gegoogelt.
Die GIS-Funktionen scheinen beim SQL-Server wohl sehr beschränkt zu sein.
Ich habe vom SQL-Server leider Null Ahnung.

Das externe Umwandeln (CSV) wäre kein Problem, aber beim/im SQL-Server kann ich leider nicht weiterhelfen.

Gruß Steffen
 
OP
H

Higgins

Geonewbie
Ich habe irgendwo im Internet eine Excel Datei gefunden, die die Umwandlung macht. Leider mit diversen Zwischenschritten, so dass es eine Heidenarbeit ist, das in SQL nachzubauen, zumal man es rückwärts machen muss, weils eben aus Excel kommt.

Sollte aber morgen fertig werden.

Ein Problem hab ich heute kurz vor Feierabend noch entdeckt, muss ich mir morgen genauer anschauen: im Hintergrund werden wohl bei den Berechnungen ein paar Nachkommastellen abgeschnitten, so dass das Ergebnis leicht abweicht. Evtl ist das am Ende aber unbedeutend genug, dass es nichts ausmacht.

Werde dann berichten.
 

steffen0815

Geocacher
Hallo,
falls du noch einen Ansatz brauchst, hier ein passender VB-Code Genauigkeit im dm-Bereich
Code:
Function LatLon(HW As Double, RW As Double, h As Single) As String
Dim pi, a, b, e, n, AL, be, ga, de, ep, e2, x, y, z, x1, y1, z1, s, t, L, B1, h1, L1, L0, B3, Bogenlänge
Dim i, BL1, BL2, BL3, BL4, Hochwert, RW1, RW2, RW3, RW4, Rechtswert, y0, b0, nf, tf, bf, pif
Dim tf1, tf2, tf3, tf4, bw, l2, l3, l4, lw, lat As Double, lon As Double
'Bessel-Ellipsoid
pi = 3.14159265358979
a = 6377397.155
b = 6356078.962
h = 4.21
e = (a ^ 2 - b ^ 2) / (a * a)
n = (a - b) / (a + b)  'nicht a*b???
AL = (a + b) / 2 * (1 + 1 / 4 * n ^ 2 + 1 / 64 * n ^ 4)
be = 3 / 2 * n - 27 / 32 * n ^ 3 + 269 / 512 * n ^ 5
ga = 21 / 16 * n ^ 2 - 55 / 32 * n ^ 4
de = 151 / 96 * n ^ 3 - 417 / 128 * n ^ 5
ep = 1097 / 512 * n ^ 4
y0 = Int(RW / 1000000)
L0 = y0 * 3
y = RW - y0 * 1000000 - 500000
b0 = HW / AL
bf = (b0 + be * Sin(2 * b0) + ga * Sin(4 * b0) + de * Sin(6 * b0) + ep * Sin(8 * b0))
nf = a / Sqr(1 - e * Sin(bf) ^ 2)
pif = Sqr(a ^ 2 / b ^ 2 * e * Cos(bf) ^ 2)
tf = Tan(bf)
tf1 = tf / 2 / nf ^ 2 * (-1 - pif ^ 2) * y ^ 2
tf2 = tf / 24 / nf ^ 4 * (5 + 3 * tf ^ 2 + 6 * pif ^ 2 - 6 * tf ^ 2 * pif ^ 2 - 4 * pif ^ 4 - 9 * tf ^ 2 * pif ^ 4) * y ^ 4
tf3 = tf / 720 / nf ^ 6 * (-61 - 90 * tf ^ 2 - 45 * tf ^ 4 - 107 * pif ^ 2 + 162 * tf ^ 2 * pif ^ 2 + 45 * tf ^ 4 * pif ^ 2) * y ^ 6
tf4 = tf / 40320 / nf ^ 8 * (1385 + 3663 * tf ^ 2 + 4095 * tf ^ 4 + 1575 * tf ^ 6) * y ^ 8
bw = (bf + tf1 + tf2) * 180 / pi
L1 = 1 / nf / Cos(bf) * y
l2 = 1 / 6 / nf ^ 3 / Cos(bf) * (-1 - 2 * tf ^ 2 - pif ^ 2) * y ^ 3
l3 = 1 / 120 / nf ^ 5 / Cos(bf) * (5 + 28 * tf ^ 2 + 24 * tf ^ 4 + 6 * pif ^ 2 + 8 * tf ^ 2 * pif ^ 2) * y ^ 5
l4 = 1 / 15040 / nf ^ 7 / Cos(bf) * (-61 - 622 * tf ^ 2 - 1320 * tf ^ 4 - 720 * tf ^ 6) * y ^ 7
lw = L0 + (L1 + l2) * 180 / pi
a = 6377397.155
b = 6356078.962
e2 = (a ^ 2 - b ^ 2) / a ^ 2
n = a / Sqr(1 - e2 * Sin(bw / 180 * pi) ^ 2)
x = (n + h) * Cos(bw / 180 * pi) * Cos(lw / 180 * pi)
y = (n + h) * Cos(bw / 180 * pi) * Sin(lw / 180 * pi)
z = (n * b ^ 2 / a ^ 2 + h) * Sin(bw / 180 * pi)
x1 = x * 1 + y * 0.0000119021759 + z * 0.000000218166156
y1 = x * -0.0000119021759 + y * 1 + z * -0.000000979323636
z1 = x * -0.000000218166156 + y * 0.0000009793236 + z * 1

x = x1 * 0.9999933 + (598.095)
y = y1 * 0.9999933 + (73.707)
z = z1 * 0.9999933 + (418.197)
a = 6378137
b = 6356752.314
e2 = (a ^ 2 - b ^ 2) / a ^ 2
s = Sqr(x ^ 2 + y ^ 2)
t = Atn(z * a / (s * b))
b = Atn((z + e2 * a ^ 2 / b * Sin(t) ^ 3) / (s - e2 * a * Cos(t) ^ 3))
L = Atn(y / x)
n = a / Sqr(1 - e2 * Sin(b) ^ 2)
h = s / Cos(b) - n
lat = b * 180 / pi
lon = L * 180 / pi
LatLon = lat & "RW" & lon

End Function

Gruß Steffen
 
OP
H

Higgins

Geonewbie
Ich habe jetzt folgendes gebaut:

Code:
 CREATE FUNCTION [dbo].[GKzuWGS84] (@Rechts FLOAT, @Hoch FLOAT) RETURNS nvarchar(255)
    BEGIN

DECLARE	@Ausgabe nvarchar(255)
		,@Hochwert FLOAT
		,@Rechtswert FLOAT
		,@transParam_dx FLOAT
		,@transParam_dy FLOAT		
		,@transParam_dz FLOAT
		,@transParam_ex FLOAT
		,@transParam_ey FLOAT
		,@transParam_ez FLOAT
		,@transParam_m FLOAT
		,@transParam_dx_dreh FLOAT
		,@transParam_dy_dreh FLOAT
		,@transParam_dz_dreh FLOAT
		,@transParam_ex_dreh FLOAT
		,@transParam_ey_dreh FLOAT
		,@transParam_ez_dreh FLOAT
		,@transParam_m_dreh FLOAT
		,@BesselEllipsoid_a FLOAT
		,@BesselEllipsoid_b FLOAT
		,@BesselEllipsoid_e FLOAT
		,@BesselEllipsoid_n FLOAT
		,@BesselEllipsoid_alpha FLOAT
		,@BesselEllipsoid_beta FLOAT
		,@BesselEllipsoid_gamma FLOAT
		,@BesselEllipsoid_delta FLOAT
		,@BesselEllipsoid_epsilon FLOAT
		,@elliptKoord_B FLOAT
		,@elliptKoord_L FLOAT
		,@elliptKoord_h FLOAT
		,@elliptKoord_N FLOAT
		,@elliptKoord_x FLOAT
		,@elliptKoord_y FLOAT
		,@elliptKoord_z FLOAT
		,@GKnBL_y0 FLOAT
		,@GKnBL_L0 FLOAT
		,@GKnBL_y FLOAT
		,@GKnBL_B0 FLOAT
		,@GKnBL_BF FLOAT
		,@GKnBL_NF FLOAT
		,@GKnBL_nanoF FLOAT
		,@GKnBL_tF FLOAT
		,@GKnBL_B1 FLOAT
		,@GKnBL_B2 FLOAT
		,@GKnBL_B3 FLOAT
		,@GKnBL_B4 FLOAT
		,@GKnBL_B FLOAT
		,@GKnBL_L1 FLOAT
		,@GKnBL_L2 FLOAT
		,@GKnBL_L3 FLOAT
		,@GKnBL_L4 FLOAT
		,@GKnBL_L FLOAT
		,@DHDNVektor1 FLOAT
		,@DHDNVektor2 FLOAT
		,@DHDNVektor3 FLOAT
		,@drehMatrixA1 FLOAT
		,@drehMatrixA2 FLOAT
		,@drehMatrixA3 FLOAT
		,@drehMatrixB1 FLOAT
		,@drehMatrixB2 FLOAT
		,@drehMatrixB3 FLOAT
		,@drehMatrixC1 FLOAT
		,@drehMatrixC2 FLOAT
		,@drehMatrixC3 FLOAT
		,@DHDNVektor1rotiert FLOAT
		,@DHDNVektor2rotiert FLOAT
		,@DHDNVektor3rotiert FLOAT
		,@DHDNVektor1massstab FLOAT
		,@DHDNVektor2massstab FLOAT
		,@DHDNVektor3massstab FLOAT
		,@DHDNVektor1translation FLOAT
		,@DHDNVektor2translation FLOAT
		,@DHDNVektor3translation FLOAT
		,@ETRF89Vektor_x FLOAT
		,@ETRF89Vektor_y FLOAT
		,@ETRF89Vektor_z FLOAT
		,@ETRF89Vektor_s FLOAT
		,@ETRF89Vektor_T FLOAT
		,@ETRF89Vektor_B FLOAT
		,@ETRF89Vektor_L FLOAT
		,@ETRF89Vektor_N FLOAT
		,@ETRF89Vektor_h FLOAT
		,@WGS84_a FLOAT
		,@WGS84_b FLOAT
		,@WGS84_e FLOAT
		,@Lat FLOAT
		,@Lon FLOAT

SELECT	@Hochwert=@Hoch
		,@Rechtswert=@Rechts
		,@transParam_dx=590.5
		,@transParam_dy=69.5
		,@transParam_dz=411.6
		,@transParam_ex=0.796
		,@transParam_ey=0.052
		,@transParam_ez= 3.601
		,@transParam_m=8.3
		,@transParam_ex_dreh=@transParam_ex*PI()/180/3600
		,@transParam_ey_dreh=@transParam_ey*PI()/180/3600
		,@transParam_ez_dreh=@transParam_ez*PI()/180/3600
		,@BesselEllipsoid_a=6377397.155
		,@BesselEllipsoid_b=6356078.962
		,@BesselEllipsoid_e=(POWER(@BesselEllipsoid_a,2)-POWER(@BesselEllipsoid_b,2))/POWER(@BesselEllipsoid_a,2)
		,@BesselEllipsoid_n=(@BesselEllipsoid_a-@BesselEllipsoid_b)/(@BesselEllipsoid_a+@BesselEllipsoid_b)
		,@BesselEllipsoid_alpha=(@BesselEllipsoid_a+@BesselEllipsoid_b)/2*(1+POWER(@BesselEllipsoid_n,2)/4+POWER(@BesselEllipsoid_n,4)/64)
		,@BesselEllipsoid_beta=(3*@BesselEllipsoid_n/2)-(27*POWER(@BesselEllipsoid_n,3)/32)+(269*POWER(@BesselEllipsoid_n,5)/512)
		,@BesselEllipsoid_gamma=(21*POWER(@BesselEllipsoid_n,2)/16)-(55*POWER(@BesselEllipsoid_n,4)/32)
		,@BesselEllipsoid_delta=(151*POWER(@BesselEllipsoid_n,3)/96)-(417*POWER(@BesselEllipsoid_n,5)/128)
		,@BesselEllipsoid_epsilon=1097*POWER(@BesselEllipsoid_n,4)/512
		,@GKnBL_y0=CAST(CAST(@Rechtswert/1000000 AS INT) AS FLOAT)
		,@GKnBL_L0=@GKnBL_y0*3
		,@GKnBL_y=@Rechtswert-@GKnBL_y0*1000000-500000
		,@GKnBL_B0=@Hochwert/@BesselEllipsoid_alpha
		,@GKnBL_BF=(@GKnBL_B0+@BesselEllipsoid_beta*SIN(2*@GKnBL_B0)+@BesselEllipsoid_gamma*SIN(4*@GKnBL_B0)+@BesselEllipsoid_delta*SIN(6*@GKnBL_B0)+@BesselEllipsoid_epsilon*SIN(8*@GKnBL_B0))
		,@GKnBL_NF=@BesselEllipsoid_a/SQRT(1-@BesselEllipsoid_e*POWER(SIN(@GKnBL_BF),2))
		,@GKnBL_nanoF=SQRT(POWER(@BesselEllipsoid_a,2)/POWER(@BesselEllipsoid_b,2)*@BesselEllipsoid_e*POWER(COS(@GKnBL_BF),2))
		,@GKnBL_tF=TAN(@GKnBL_BF)
		,@GKnBL_B1=@GKnBL_tF/2/POWER(@GKnBL_NF,2)*(-1-POWER(@GKnBL_nanoF,2))*POWER(@GKnBL_y,2)
		,@GKnBL_B2=@GKnBL_tF/24/POWER(@GKnBL_NF,4)*(5+3*POWER(@GKnBL_tF,2)+6*POWER(@GKnBL_nanoF,2)-6*POWER(@GKnBL_tF,2)*POWER(@GKnBL_nanoF,2)-4*POWER(@GKnBL_nanoF,4)-9*POWER(@GKnBL_tF,2)*POWER(@GKnBL_nanoF,4))*POWER(@GKnBL_y,4)
		,@GKnBL_B3=@GKnBL_tF/720/POWER(@GKnBL_NF,6)*(-61-90*POWER(@GKnBL_tF,2)-45*POWER(@GKnBL_tF,4)-107*POWER(@GKnBL_nanoF,2)+162*POWER(@GKnBL_tF,2)*POWER(@GKnBL_nanoF,2)+45*POWER(@GKnBL_tF,4)*POWER(@GKnBL_nanoF,2))*POWER(@GKnBL_y,6)
		,@GKnBL_B4=@GKnBL_tF/40320/POWER(@GKnBL_NF,8)*(1385+3663*POWER(@GKnBL_tF,2)+4095*POWER(@GKnBL_tF,4)+1575*POWER(@GKnBL_tF,6))*POWER(@GKnBL_y,8)
		,@GKnBL_B=(@GKnBL_BF+@GKnBL_B1+@GKnBL_B2)*180/PI()
		,@GKnBL_L1=1/@GKnBL_NF/COS(@GKnBL_BF)*@GKnBL_y
		,@GKnBL_L2=1/POWER(@GKnBL_NF,3)/6/COS(@GKnBL_BF)*(-1-2*POWER(@GKnBL_tF,2)-POWER(@GKnBL_nanoF,2))*POWER(@GKnBL_y,3)
		,@GKnBL_L3=1/POWER(@GKnBL_NF,5)/120/COS(@GKnBL_BF)*(5+28*POWER(@GKnBL_tF,2)+24*POWER(@GKnBL_tF,4)+6*POWER(@GKnBL_nanoF,2)+8*POWER(@GKnBL_tF,2)*POWER(@GKnBL_nanoF,2))*POWER(@GKnBL_y,5)
		,@GKnBL_L4=1/POWER(@GKnBL_NF,7)/15040/COS(@GKnBL_BF)*(-61-622*POWER(@GKnBL_tF,2)-1320*POWER(@GKnBL_tF,4)-720*POWER(@GKnBL_tF,6))*POWER(@GKnBL_y,7)
		,@GKnBL_L=@GKnBL_L0+((@GKnBL_L1+@GKnBL_L2)*180/PI())
		,@elliptKoord_B=@GKnBL_B
		,@elliptKoord_L=@GKnBL_L
		,@elliptKoord_h=6.0
		,@elliptKoord_N=@BesselEllipsoid_a/SQRT(1-(@BesselEllipsoid_e*(SIN(@elliptKoord_B/180*PI())*SIN(@elliptKoord_B/180*PI()))))
		,@elliptKoord_x=(@elliptKoord_N+@elliptKoord_h)*COS(@elliptKoord_B/180*PI())*COS(@elliptKoord_L/180*PI())
		,@elliptKoord_y=(@elliptKoord_N+@elliptKoord_h)*COS(@elliptKoord_B/180*PI())*SIN(@elliptKoord_L/180*PI())
		,@elliptKoord_z=((@elliptKoord_N*POWER(@BesselEllipsoid_b,2)/POWER(@BesselEllipsoid_a,2))+@elliptKoord_h)*SIN(@elliptKoord_B/180*PI())
		,@drehMatrixA1=1
		,@drehMatrixB1=@transParam_ez_dreh
		,@drehMatrixC1=@transParam_ey_dreh
		,@drehMatrixA2=@drehMatrixB1*(-1)
		,@drehMatrixB2=1	
		,@drehMatrixC2=@transParam_ex_dreh
		,@drehMatrixA3=@drehMatrixC1*(-1)
		,@drehMatrixB3=@drehMatrixC2*(-1)
		,@drehMatrixC3=1
		,@transParam_m_dreh=1-(@transParam_m*0.000001)
		,@transParam_dx_dreh=((@transParam_dx*@drehMatrixA1)+(@transParam_dy*@drehMatrixB1)+(@transParam_dz*@drehMatrixC1))*@transParam_m_dreh
		,@transParam_dy_dreh=((@transParam_dx*@drehMatrixA2)+(@transParam_dy*@drehMatrixB2)+(@transParam_dz*@drehMatrixC2))*@transParam_m_dreh
		,@transParam_dz_dreh=((@transParam_dx*@drehMatrixA3)+(@transParam_dy*@drehMatrixB3)+(@transParam_dz*@drehMatrixC3))*@transParam_m_dreh		
		,@DHDNVektor1=@elliptKoord_x
		,@DHDNVektor2=@elliptKoord_y
		,@DHDNVektor3=@elliptKoord_z
		,@DHDNVektor1rotiert=(@DHDNVektor1*@drehMatrixA1)+(@DHDNVektor2*@drehMatrixB1)+(@DHDNVektor3*@drehMatrixC1)
		,@DHDNVektor2rotiert=(@DHDNVektor1*@drehMatrixA2)+(@DHDNVektor2*@drehMatrixB2)+(@DHDNVektor3*@drehMatrixC2)
		,@DHDNVektor3rotiert=(@DHDNVektor1*@drehMatrixA3)+(@DHDNVektor2*@drehMatrixB3)+(@DHDNVektor3*@drehMatrixC3)
		,@DHDNVektor1massstab=@DHDNVektor1rotiert*@transParam_m_dreh
		,@DHDNVektor2massstab=@DHDNVektor2rotiert*@transParam_m_dreh
		,@DHDNVektor3massstab=@DHDNVektor3rotiert*@transParam_m_dreh
		,@DHDNVektor1translation=@DHDNVektor1massstab+@transParam_dx_dreh
		,@DHDNVektor2translation=@DHDNVektor2massstab+@transParam_dy_dreh
		,@DHDNVektor3translation=@DHDNVektor3massstab+@transParam_dz_dreh
		,@WGS84_a=6378137.000
		,@WGS84_b=6356752.314
		,@WGS84_e=(POWER(@WGS84_a,2)-POWER(@WGS84_b,2))/POWER(@WGS84_a,2)
		,@ETRF89Vektor_x=@DHDNVektor1translation
		,@ETRF89Vektor_y=@DHDNVektor2translation
		,@ETRF89Vektor_z=@DHDNVektor3translation
		,@ETRF89Vektor_s=SQRT(POWER(@ETRF89Vektor_x,2)+POWER(@ETRF89Vektor_y,2))
		,@ETRF89Vektor_T=ATAN(@ETRF89Vektor_z*@WGS84_a/(@ETRF89Vektor_s*@WGS84_b))
		,@ETRF89Vektor_B=ATAN((@ETRF89Vektor_z+@WGS84_e*@WGS84_a*@WGS84_a/@WGS84_b*POWER(SIN(@ETRF89Vektor_T),3))/(@ETRF89Vektor_s-@WGS84_e*@WGS84_a*POWER(COS(@ETRF89Vektor_T),3)))
		,@ETRF89Vektor_L=ATAN(@ETRF89Vektor_y/@ETRF89Vektor_x)
		,@ETRF89Vektor_N=@WGS84_a/SQRT(1-@WGS84_e*POWER(SIN(@ETRF89Vektor_B),2))
		,@ETRF89Vektor_h=@ETRF89Vektor_s/COS(@ETRF89Vektor_B)-@ETRF89Vektor_N
		,@Lat=@ETRF89Vektor_B*180/PI()
		,@Lon=@ETRF89Vektor_L*180/PI()
		,@Ausgabe=CONCAT(CAST(CAST(@Lat AS decimal(25,10)) AS nvarchar(100)),' ',CAST(CAST(@Lon AS decimal(25,10)) AS nvarchar(100)))

RETURN @Ausgabe
END

GO

Rückwärts aus der Excel-Tabelle abgeleitet, die ich gefunden hatte. Funktioniert soweit gut, ein bisschen ungenau in den Nachkommastellen, weil SQL bei den Zwischenrechnungen ein paar abschneidet, wirkt sich im Ergebnis aber um maximal 3m oder so aus, was für meine Zwecke genau genug ist. h habe ich von 4.21 auf 6 gesetzt, weil meine Orte alle in Hamburg liegen und Hamburg laut Wikipedia durchschnittlich 6m über Normalnull liegt...
 
Oben