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?
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
No comments:
Post a Comment