0

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?

  • Are you wanting these added to new fields? Do you have any GIS software - potentially you could use that to do the conversion for you – Midavalo Oct 20 '16 at 4:29
  • I don't use SQL Server, but is there a method to do the coordinate transformation in the database as part of a query? E.g. with PostGIS I'd use ST_Transform. See gis.stackexchange.com/questions/19917/… – Richard Law Oct 20 '16 at 4:38
  • @midavalo we do have gis but i want a fast and dynamic way to convert this. The data is mastered in a non GIS application and it does shift. As a rule everything we do is in NZTM. Im only looking at this for compatibility with bing in Excel as a proof of concept for an adhoc analysis tool. – Sir Swears-a-lot Oct 20 '16 at 7:48
1

The following SQL code is due to Peter Moore at Gisborne District Council. I had merely provided the conversion formulae in Excel following the algorithms described by LINZ.

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

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.