[GIS] Convert Latitude/Longitude to Lambert Conformal Conic in SQL

coordinate systemlambert-conformal-conicsql

How can I convert the from Latitude/Longitude to Lambert Conic Conformal using SQL?

I've been working on project that will require me to convert between NZTM And LCC projections. It's related to this post.

LINZ don't provide a formula to convert directly between NZTM and LCC, but both can be converted to/from Lat/Long. I've been very fortunate in getting some help with the NZTM conversion and am now working on the LCC conversion.

Best Answer

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
Related Question