The following SQL code is an implementation of the projection conversion formulas provided by LINZ. I'm sharing here in the hope it may help someone else.
To convert from Lat/Long to LCC:
Begin
Declare
-- constants
@GRS80 float = 6378137,
@InverseFlattening float = 298.2572221,
@Firststandardparallel float = -37.5 ,
@Secondstandardparallel float = -44.5 ,
@Originlatitude float = -41.0 ,
@Originlongitude float = 173.0 ,
@FalseNorthing float = 7000000.0 ,
@FalseEasting float = 3000000.0 ,
--input Latitude (Y), Longitude (X)
@Y float = 41.72,
@X float = 172.48,
@a float = 0,
@f float = 0,
@phi1 float = 0,
@phi2 float = 0,
@phi0 float = 0,
@lambda0 float = 0,
@N0 float = 0,
@E0 float = 0,
@phiIn float = 0,
@lamdaIn float = 0,
@e float = 0,
@m float = 0,
@t float = 0,
@m1 float = 0,
@m2 float = 0,
@t0 float = 0,
@t1 float = 0,
@t2 float = 0,
@n float = 0,
@Fcap float = 0,
@rho float = 0,
@rho0 float = 0,
@gamma float = 0,
@Nout float = 0,
@Eout float = 0;
set @a = @GRS80
set @f = 1/@InverseFlattening
set @phi1 = @Firststandardparallel * pi()/180
set @phi2 = @Secondstandardparallel * pi()/180
set @phi0 = @Originlatitude * pi()/180
set @lambda0 = @Originlongitude * pi()/180
set @N0 = @FalseNorthing
set @E0 = @FalseEasting
set @phiIn = -@Y * pi()/180
set @lamdaIn = @X * pi()/180
set @e = sqrt(2*@f-power(@f,2))
set @m = cos(@phiIn)/sqrt(1-power((@e*SIN(@phiIn)),2))
set @t = tan(0.25*pi()-.5*@phiIn)*power((1+@e*sin(@phiIn)),.5*@e)/power((1-@e*sin(@phiIn)),0.5*@e)
set @m1 = cos(@phi1)/sqrt(1-power(@e*SIN(@phi1),2))
set @m2 = cos(@phi2)/sqrt(1-power(@e*SIN(@phi2),2))
set @t0 = tan(pi()/4-@phi0/2)/power((1-@e*SIN(@phi0))/(1+@e*SIN(@phi0)),(@e/2))
set @t1 = tan(pi()/4-@phi1/2)/power((1-@e*SIN(@phi1))/(1+@e*SIN(@phi1)),(@e/2))
set @t2 = tan(pi()/4-@phi2/2)/power((1-@e*SIN(@phi2))/(1+@e*SIN(@phi2)),(@e/2))
set @n = (LOG(@m1)-LOG(@m2))/(LOG(@t1)-LOG(@t2))
set @Fcap = @m1/(@n*power(@t1,@n))
set @rho = @a*@Fcap*power(@t,@n)
set @rho0 = @a*@Fcap*power(@t0,@n)
set @gamma = @n*(@lamdaIn-@lambda0)
set @Nout = @N0+@rho0-@rho*cos(@gamma)
set @Eout = @E0+@rho*SIN(@gamma) ;
print N' @Nout : ' + str( @Nout ,10,8)
print N' @Eout : ' + str( @Eout ,10,8)
End
To Convert back:
Begin
Declare
-- constants
@GRS80 float = 6378137,
@InverseFlattening float = 298.2572221,
@Firststandardparallel float = -37.5 ,
@Secondstandardparallel float = -44.5 ,
@Originlatitude float = -41.0 ,
@Originlongitude float = 173.0 ,
@FalseNorthing float = 7000000.0 ,
@FalseEasting float = 3000000.0 ,
--input Northing (Y), Easting (X)
@Nin float = 6920054.29,
@Ein float = 2956806.82,
-- variables.
@a float = 0,
@f float = 0,
@phi1 float = 0,
@phi2 float = 0,
@phi0 float = 0,
@lambda0 float = 0,
@N0 float = 0,
@E0 float = 0,
@lamdaIn float = 0,
@e float = 0,
@m float = 0,
@t float = 0,
@m1 float = 0,
@m2 float = 0,
@t0 float = 0,
@t1 float = 0,
@t2 float = 0,
@n float = 0,
@Fcap float = 0,
@rho float = 0,
@rho0 float = 0,
@gamma float = 0,
@Nprime float = 0,
@Eprime float = 0,
@rhoprime float = 0,
@tprime float = 0,
@gammaprime float = 0,
@phiout float = 0 ,
@phix float = 0 ,
@LatOut float = 0 ,
@LongOut float = 0
;
SET @a = @GRS80
SET @f = 1/@InverseFlattening
SET @phi1 = @Firststandardparallel * PI()/180
SET @phi2 = @Secondstandardparallel * PI()/180
SET @phi0 = @Originlatitude * PI()/180
SET @lambda0 = @Originlongitude * PI()/180
SET @N0 = @FalseNorthing
SET @E0 = @FalseEasting
SET @e = sqrt(2*@f-power(@f,2))
SET @m1 = COS(@phi1)/SQRT(1-power(@e*SIN(@phi1),2))
SET @m2 = COS(@phi2)/SQRT(1-power(@e*SIN(@phi2),2))
SET @t0 = TAN(PI()/4-@phi0/2)/power((1-@e*SIN(@phi0))/(1+@e*SIN(@phi0)),(@e/2))
SET @t1 = TAN(PI()/4-@phi1/2)/power((1-@e*SIN(@phi1))/(1+@e*SIN(@phi1)),(@e/2))
SET @t2 = TAN(PI()/4-@phi2/2)/power((1-@e*SIN(@phi2))/(1+@e*SIN(@phi2)),(@e/2))
SET @n = (LOG(@m1)-LOG(@m2))/(LOG(@t1)-LOG(@t2))
SET @Fcap = @m1/(@n*power(@t1,@n))
SET @rho0 = @a*@Fcap*power(@t0,@n)
SET @gamma = @n*(@lamdaIn-@lambda0)
SET @Nprime = @Nin-@N0
SET @Eprime = @Ein-@E0
SET @rhoprime = SIGN(@n)*SQRT(power(@Eprime,2)+power((@rho0-@Nprime),2))
SET @tprime = power((@rhoprime/(@a*@Fcap)),(1/@n))
SET @gammaprime = ATAN(@Eprime/(@rho0-@Nprime))
SET @phiout = PI()/2-2*ATAN(@tprime)
SET @phix = dbo.fn_phiLoop(@phiout, @tprime , @e)
SET @LatOut = @phix *180/PI()
SET @LongOut = (@gammaprime/@n+@lambda0)*180/PI() ;
print N' @Latitude : ' + STR( @LatOut ,10,8)
print N' @Longitude : ' + STR( @LongOut ,10,8)
End
The 2nd formula uses an iterative loop to achieve convergence.
CREATE FUNCTION [dbo].[fn_phiLoop](@phiIN float, @tprime float, @e float )
RETURNS float
AS
BEGIN
DECLARE
@phiOUT float = 0,
@cnt int =0;
WHILE (@cnt < 10)
BEGIN
SET @phiOUT = PI()/2-2*ATAN( @tprime*power(((1-@e*SIN(@phiIN))/((1+@e*SIN(@phiIN)))) , (@e/2)))
SET @phiIN = @phiOUT
SET @cnt = @cnt + 1
END
RETURN @phiOUT ;
END
GO