0

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.

|improve this question
  • SQL Server doesn't support coordinate system transformations. There are other questions and hacks outlined on this forum. – DPSSpatial Mar 26 '18 at 1:49
  • I'm a DBA. When the only tool you've got is a hammer... everything looks like a nail ;) – Sir Swears-a-lot Mar 26 '18 at 2:50
  • I hear ya... we use MS SQL Server, but maintain then ETL our data in PostGIS to MSSQL... works like a charm, but not sure if you have that need... – DPSSpatial Mar 26 '18 at 14:30
  • @DPSSpatial we run ArcGIS with ArcSDE on SQL. But we do everything in NZTM. I am working with NetCDF files based in LCC and need to translate points between them. – Sir Swears-a-lot Mar 27 '18 at 7:27
1

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
|improve this answer

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.