SQL Server – Converting NZTM to Latitude/Longitude

coordinate systemlatitude longitudesql server

I have NZTM coordinate data stored in SQL Server in a Geometry field.

I'd like to convert these coordinates to Lat/Long so that it can be used with Bing maps in Excel (and a few other things).

I have found the formula to convert from NZTM to Lat/Long here on the LINZ website.

I'm not a mathematician, and that formula looks horrific.

Can anyone share a working example of the code in VBA or SQL that I could use to convert this on the fly?

Best Answer

The following SQL query is an implementation of the formula provided by LINZ to convert from NZTM to Geographic (Lat/Long). Source

This example uses the centre of the amphitheater at Te Papa as a reference point.

Input: Northing 5427502.0, Easting 1749165.0 Output: Longitude -41.290153, Latitude 174.781428

Begin 

Declare
-- Constants
@a      float = 6378137,
@f      float = 1/298.257222101,
@phizero float =    0,
@lambdazero float = 173,
@Nzero  float = 10000000,
@Ezero  float = 1600000,
@kzero  float = 0.9996,

-- Common variables
@b  float = 0,
@esq float  = 0,
@mzero float = 0,
@A0 float = 0,
@A2 float = 0,
@A4 float = 0,
@A6 float = 0,

--input Northing (Y), Easting (X) variables 
@N      float = 5427502.0, 
@E      float = 1749165.0, 

-- Variables
@Nprime float = 0,
@mprime float = 0,
@smn    float = 0,
@G      float = 0,
@sigma  float = 0,
@phiprime   float = 0,
@rhoprime   float = 0,
@upsilonprime float = 0,
@psiprime   float = 0,
@tprime float = 0,
@Eprime float = 0,
@chi    float = 0,
@term_1 float = 0,
@term_2 float = 0,
@term_3 float = 0,
@term_4 float = 0,
@term1  float = 0,
@term2  float = 0,
@term3  float = 0,
@term4  float = 0,

--Output
@latitude   float = 0, 
@longitude  float = 0; 

--Calculation: To NZTM to NZGD2000/WGS84

--Common

SET @b          = @a*(1-@f)
SET @esq        = 2*@f-POWER(@f,2)
SET @A0         =1-@esq/4-3*POWER(@esq,2)/64-5*POWER(@esq,3)/256 
SET @A2         =0.375*(@esq+POWER(@esq,2)/4+15*power(@esq,3)/128) 
SET @A4         =15*(POWER(@esq,2)+3*POWER(@esq,3)/4)/256  
SET @A6         =35*POWER(@esq,3)/3072 

SET @Nprime     = @N-@Nzero
SET @mprime     = @Nprime/@kzero
SET @smn        = (@a-@b)/(@a+@b)
SET @G          = @a*(1-@smn)*(1-POWER(@smn,2))*(1+9*POWER(@smn,2)/4+225*POWER(@smn,4)/64)*PI()/180.0
SET @sigma      = @mprime*PI()/(180*@G)
SET @phiprime   = @sigma+(3*@smn/2-27*POWER(@smn,3)/32)*SIN(2*@sigma)+(21*POWER(@smn,2)/16-55*POWER(@smn,4)/32)*SIN(4*@sigma)+(151*POWER(@smn,3)/96)*SIN(6*@sigma)+(1097*POWER(@smn,4)/512)*SIN(8*@sigma)
SET @rhoprime   = @a*(1-@esq)/POWER((1-@esq*POWER((SIN(@phiprime)),2)),1.5)
SET @upsilonprime   =@a/SQRT(1-@esq*POWER((SIN(@phiprime)),2)) 

SET @psiprime   = @upsilonprime/@rhoprime
SET @tprime     = TAN(@phiprime)
SET @Eprime     = @E-@Ezero
SET @chi        = @Eprime/(@kzero*@upsilonprime)
SET @term_1     = @tprime*@Eprime*@chi/(@kzero*@rhoprime*2)
SET @term_2     = @term_1*POWER(@chi,2)/12*(-4*POWER(@psiprime,2)+9*@psiprime*(1-POWER(@tprime,2))+12*POWER(@tprime,2))
SET @term_3     = @tprime*@Eprime*POWER(@chi,5)/(@kzero*@rhoprime*720)*(8*POWER(@psiprime,4)*(11-24*POWER(@tprime,2))-12*POWER(@psiprime,3)*(21-71*POWER(@tprime,2))+15*POWER(@psiprime,2)*(15-98*POWER(@tprime,2)+15*POWER(@tprime,4))+180*@psiprime*(5*POWER(@tprime,2)-3*POWER(@tprime,4))+360*POWER(@tprime,4))
SET @term_4     = @tprime*@Eprime*POWER(@chi,7)/(@kzero*@rhoprime*40320)*(1385+3633*POWER(@tprime,2)+4095*POWER(@tprime,4)+1575*POWER(@tprime,6))
SET @term1      = @chi*(1/COS(@phiprime))
SET @term2      = POWER(@chi,3)*(1/COS(@phiprime))/6*(@psiprime+2*POWER(@tprime,2))
SET @term3      = POWER(@chi,5)*(1/COS(@phiprime))/120*(-4*POWER(@psiprime,3)*(1-6*POWER(@tprime,2))+POWER(@psiprime,2)*(9-68*POWER(@tprime,2))+72*@psiprime*POWER(@tprime,2)+24*POWER(@tprime,4))
SET @term4      = POWER(@chi,7)*(1/COS(@phiprime))/5040*(61+662*POWER(@tprime,22)+1320*POWER(@tprime,4)+720*POWER(@tprime,6))

SET @latitude   = (@phiprime-@term_1+@term_2-@term_3+@term_4)*180/PI()
SET @longitude  = @lambdazero+180/PI()*(@term1-@term2+@term3-@term4)

print N' @a:        ' + STR( @a         ,10,8) 
print N' @b:        ' + STR( @b         ,10,8) 
print N' @esq:      ' + STR( @esq       ,10,8) 
print N' @A0:       ' + STR( @A0        ,10,8) 
print N' @A2:       ' + STR( @A2        ,10,8) 
print N' @A4:       ' + STR( @A4        ,10,8) 
print N' @A6:       ' + STR( @A6        ,15,13) 
print N' @Nprime:   ' + STR( @Nprime    ,10,8) 
print N' @mprime:   ' + STR( @mprime    ,10,8) 
print N' @smn:      ' + STR( @smn       ,10,8) 
print N' @G:        ' + STR( @G         ,10,8) 
print N' @sigma:    ' + STR( @sigma     ,10,8) 
print N' @phiprime: ' + STR( @phiprime  ,10,8) 
print N' @rhoprime: ' + STR( @rhoprime  ,10,8) 
print N' @upsilonprime: ' + STR( @upsilonprime ,10,8) 
print N' @psiprime: ' + STR( @psiprime  ,10,8) 
print N' @tprime:   ' + STR( @tprime    ,10,8) 
print N' @Eprime:   ' + STR( @Eprime    ,10,8) 
print N' @chi:      ' + STR( @chi       ,10,8) 
print N' @term_1:   ' + STR( @term_1    ,10,8) 
print N' @term_2:   ' + STR( @term_2    ,10,8) 
print N' @term_3:   ' + STR( @term_3    ,10,8) 
print N' @term_4:   ' + STR( @term_4    ,10,8) 
print N' @term1:    ' + STR( @term1     ,10,8) 
print N' @term2:    ' + STR( @term2     ,10,8) 
print N' @term3:    ' + STR( @term3     ,10,8) 
print N' @term4:    ' + STR( @term4     ,10,8) 

-- annnddd here we are :)   
print N' @latitude: ' + STR( @latitude  ,10,8) 
print N' @longitude:' + STR( @longitude ,10,8) 

End

Another example uses the centre of the amphitheatre at Te Papa as a reference point.

Input: Longitude -41.290153, Latitude 174.781428 Output: Northing 5427502.0, Easting 1749165.0

Begin 
Declare
-- constants
@a      float = 6378137,
@f      float = 1/298.257222101,
@phizero float =    0,
@lambdazero float = 173,
@Nzero  float = 10000000,
@Ezero  float = 1600000,
@kzero  float = 0.9996,

-- Common 
@b  float = 0,
@esq float  = 0,
@mzero float = 0,
@A0 float = 0,
@A2 float = 0,
@A4 float = 0,
@A6 float = 0,

--input Latitude (Y), Longitude (X) and varibles 

@phi    float = -41.290153, 
@lambda float = 174.781428,

-- Variables
@rho    float = 0,
@upsilon float = 0,
@psi    float = 0,
@t      float = 0,
@omega  float = 0,
@m      float = 0,
@t_1    float = 0,
@t_2    float = 0,
@t_3    float = 0,
@t_4    float = 0,
@t1     float = 0,
@t2     float = 0,
@t3     float = 0, 

--Output
@Easting    float = 0, 
@Northing   float = 0, 
@latitude   float = 0, 
@longitude  float = 0
; 

--Calculation: NZGD2000/WGS84 to NZTM

--Common

SET @b          = @a*(1-@f)
SET @esq        = 2*@f-POWER(@f,2)
SET @A0         = 1-@esq/4-3*POWER(@esq,2)/64-5*POWER(@esq,3)/256 
SET @A2         = 0.375*(@esq+POWER(@esq,2)/4+15*power(@esq,3)/128)  
SET @A4         = 15*(POWER(@esq,2)+3*POWER(@esq,3)/4)/256  
SET @A6         = 35*POWER(@esq,3)/3072     

-- Calculation  To NZTM
SET @rho        = @a*(1-@esq)/POWER(1-@esq*POWER((SIN(PI()*@phi/180)),2),1.5) 
SET @upsilon    = @a/SQRT(1-@esq*POWER((SIN(PI()*@phi/180)),2))
SET @psi        = @upsilon/@rho
SET @t          = TAN(@phi*PI()/180)
SET @omega      = (@lambda-@lambdazero)*PI()/180
SET @m          = @a*(@A0*PI()*@phi/180-@A2*SIN(2*PI()*@phi/180)+@A4*SIN(4*PI()*@phi/180)-@A6*SIN(6*PI()*@phi/180))

SET @t_1        = POWER(@omega,2)*@upsilon/2*SIN(@phi*PI()/180)*COS(@phi*PI()/180)
SET @t_2        = POWER(@omega,4)*@upsilon/24*SIN(@phi*PI()/180)*POWER((COS(@phi*PI()/180)),3)*(4*POWER(@psi,2)+@psi-POWER(@t,2))
SET @t_3        = POWER(@omega,6)*@upsilon/720*SIN(PI()*@phi/180)*POWER((COS(PI()*@phi/180)),5)*(8*POWER(@psi,4)*(11-24*POWER(@t,2))-28*POWER(@psi,3)*(1-6*POWER(@t,2))+POWER(@psi,2)*(1-32*POWER(@t,2))-@psi*2*POWER(@t,2)+POWER(@t,4))
SET @t_4        = POWER(@omega,8)*@upsilon/40320*SIN(PI()*@phi/180)*POWER((COS(PI()*@phi/180)),7)*(1385-3111*POWER(@t,2)+543*POWER(@t,4)-POWER(@t,6))
SET @t1         = POWER(@omega,2)/6*POWER((COS(PI()*@phi/180)),2)*(@psi-POWER(@t,2))
SET @t2         = POWER(@omega,4)/120*POWER((COS(PI()*@phi/180)),4)*(4*POWER(@psi,3)*(1-6*POWER(@t,2))+POWER(@psi,2)*(1+8*POWER(@t,2))-2*@psi*POWER(@t,2)+POWER(@t,4))
SET @t3         = POWER(@omega,6)/5040*POWER((COS(PI()*@phi/180)),6)*(61-479*POWER(@t,2)+179*POWER(@t,4)-POWER(@t,6))

SET @Easting    = @Ezero+@kzero*@upsilon*@omega*COS(PI()*@phi/180)*(1+@t1+@t2+@t3)
SET @Northing   = @Nzero+@kzero*(@m+@t_1+@t_2+@t_3+@t_4)

print N' @b:        ' + STR( @b         ,12,10) 
print N' @esq:      ' + STR( @esq       ,12,10) 
print N' @A0:       ' + STR( @A0        ,12,10) 
print N' @A2:       ' + STR( @A2        ,12,10) 
print N' @A4:       ' + STR( @A4        ,12,10) 
print N' @A6:       ' + STR( @A6        ,12,10) 
print N' @rho:      ' + STR( @rho       ,12,10) 
print N' @upsilon:  ' + STR( @upsilon   ,12,10) 
print N' @psi:      ' + STR( @psi       ,12,10) 
print N' @t:        ' + STR( @t         ,12,10) 
print N' @omega:    ' + STR( @omega     ,12,10) 
print N' @m:        ' + STR( @m         ,12,10) 
print N' @t_1:      ' + STR( @t_1       ,12,10) 
print N' @t_2:      ' + STR( @t_2       ,12,10) 
print N' @t_3:      ' + STR( @t_3       ,12,10) 
print N' @t_4:      ' + STR( @t_4       ,12,10) 
print N' @t1:       ' + STR( @t1        ,12,10) 
print N' @t2:       ' + STR( @t2        ,12,10) 
print N' @t3:       ' + STR( @t3        ,12,10) 

print N' @Easting:  ' + STR( @Easting   ,12,10) 
print N' @Northing: ' + STR( @Northing  ,12,10) 

End

Best Answer

Related Solutions

R Map Library – Converting NZMG or NZTM to Latitude/Longitude

[GIS] How to convert 6 & 7 digit coordinates to degrees latitude & longitude

Related Question