Skip to content

Lambert Conformal Conic (tangent case) support #224

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Lambert Conformal Conic (tangent case) support #224

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
cd279fc
Adding support to Lambert Conic Conformal projections Tangent case.
ramespada May 15, 2023
86ea456
Minor bugs corrections: I forgot to add the xyll_lam.f90 file and pas…
ramespada May 16, 2023
bdb71cf
Merge branch 'USEPA:main' into new_branch
ramespada Oct 19, 2023
45c2196
Merge branch 'USEPA:main' into proj_lcc_tan
ramespada Dec 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions PREP/mcip/src/ll2xy_lam.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,9 @@ SUBROUTINE ll2xy_lam (phi, lambda, phi1, phi2, lambda0, phi0, xx, yy)
!-------------------------------------------------------------------------------

IF ( ABS( phi1 - phi2 ) < phitol ) THEN ! tangent case
WRITE (*,f9000) TRIM(pname), phi1, phi2
CALL graceful_stop (pname)
! CALL ll2xy_lam_tan (phi, lambda, phi1, lambda0, xx, yy)
! WRITE (*,f9000) TRIM(pname), phi1, phi2
! CALL graceful_stop (pname)
CALL ll2xy_lam_tan (phi, lambda, phi1, phi2, lambda0, phi0, xx, yy)
ELSE ! secant case
CALL ll2xy_lam_sec (phi, lambda, phi1, phi2, lambda0, phi0, xx, yy)
ENDIF
Expand Down
41 changes: 28 additions & 13 deletions PREP/mcip/src/ll2xy_lam_tan.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
! subject to their copyright restrictions. !
!------------------------------------------------------------------------------!

SUBROUTINE ll2xy_lam_tan (phi, lambda, phi0, lambda0, xx, yy)
SUBROUTINE ll2xy_lam_tan (phi, lambda, phi1, phi2, lambda0, phi0, xx, yy)

!-------------------------------------------------------------------------------
! Name: Latitude-Longitude to (X,Y) for Lambert Conformal Projection
Expand All @@ -41,14 +41,22 @@ SUBROUTINE ll2xy_lam_tan (phi, lambda, phi0, lambda0, xx, yy)
REAL, INTENT(IN) :: lambda0 ! standard longitude [deg]
REAL, INTENT(IN) :: phi ! latitude [deg]
REAL(8) :: phirad ! latitude [rad]
REAL, INTENT(IN) :: phi0 ! true latitude [deg]
REAL(8) :: phi0rad ! true latitude [rad]
REAL, INTENT(IN) :: phi0 ! reference latitude [deg]
REAL(8) :: phi0rad ! reference latitude [rad]
REAL, INTENT(IN) :: phi1 ! true latitude [deg]
REAL(8) :: phi1rad ! true latitude [rad]
REAL, INTENT(IN) :: phi2 ! true latitude [deg]
!phi2 (not used)
!phi2rad (not used)
REAL(8) :: pi
REAL(8) :: piover4 ! pi/4
REAL(8) :: psi !auxiliary function
REAL(8) :: rho ! polar radius to origin
REAL(8) :: rho0 ! polar radius to latitude phi
REAL(8) :: term
REAL(8) :: term0
REAL(8) :: term1
!term2 (not used)
REAL(8) :: theta ! polar angle
REAL(8) :: sinphi0 ! cone constant
REAL, INTENT(OUT) :: xx ! X-coordinate from origin
Expand All @@ -68,36 +76,43 @@ SUBROUTINE ll2xy_lam_tan (phi, lambda, phi0, lambda0, xx, yy)
! Compute cone constant, SINPHI0.
!-------------------------------------------------------------------------------

phi0rad = phi0 * deg2rad ! convert PHI0 from degrees to radians
sinphi0 = DSIN (phi0rad)
phi0rad = DBLE(phi0) * deg2rad ! convert PHI0 from degrees to radians
phi1rad = DBLE(phi1) * deg2rad ! convert PHI1 from degrees to radians

term0 = DTAN (piover4 - phi0rad/2.0d0)
term1 = DTAN (piover4 - phi1rad/2.0d0)
!term2 (not used)


sinphi0 = DSIN (phi1rad) !this is the only difference with the secant case.

!-------------------------------------------------------------------------------
! Compute polar angle, THETA.
!-------------------------------------------------------------------------------

dlambda = (lambda - lambda0) * deg2rad
dlambda = DBLE(lambda - lambda0) * deg2rad
theta = dlambda * sinphi0

!-------------------------------------------------------------------------------
! Compute polar radius to origin, RHO0, where origin is at PHI0.
!-------------------------------------------------------------------------------

rho0 = drearth * DCOS(phi0rad) / sinphi0
psi = drearth * DCOS(phi1rad) / sinphi0 / (term1**sinphi0)
rho0 = psi * (term0**sinphi0)

!-------------------------------------------------------------------------------
! Compute polar radius to latitude PHI, RHO.
!-------------------------------------------------------------------------------

phirad = phi * deg2rad ! convert PHI from degrees to radians
term = DTAN (piover4 - phirad /2.0d0)
term0 = DTAN (piover4 - phi0rad/2.0d0)
rho = rho0 * (( term / term0 )**sinphi0)
phirad = DBLE(phi) * deg2rad ! convert PHI from degrees to radians
term = DTAN (piover4 - phirad/2.0d0)
rho = psi * (term**sinphi0)

!-------------------------------------------------------------------------------
! Compute Cartesian coordinates, XX and YY.
!-------------------------------------------------------------------------------

xx = REAL( rho * DSIN(theta))
yy = REAL(rho0 - rho * DCOS(theta))
xx = REAL( rho * DSIN(theta) )
yy = REAL( rho0 - rho * DCOS(theta) )

END SUBROUTINE ll2xy_lam_tan
10 changes: 8 additions & 2 deletions PREP/mcip/src/mapfac_lam.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,14 @@ REAL FUNCTION mapfac_lam (phi, phi1, phi2)
term1 = DTAN (piover4 - phi1rad/2.0d0)
term2 = DTAN (piover4 - phi2rad/2.0d0)

sinphi0 = DLOG ( DCOS(phi1rad) / DCOS(phi2rad) )
sinphi0 = sinphi0 / DLOG (term1 / term2)
if ( ABS(phi1 - phi2) > 0.1) then
sinphi0 = DLOG ( DCOS(phi1rad) / DCOS(phi2rad) )
sinphi0 = sinphi0 / DLOG (term2 / term1)
else
sinphi0 = DSIN(phi1rad)
endif
!sinphi0 = DLOG ( DCOS(phi1rad) / DCOS(phi2rad) )
!sinphi0 = sinphi0 / DLOG (term1 / term2)

!-------------------------------------------------------------------------------
! Compute map-scale factor, MAPFAC.
Expand Down
8 changes: 6 additions & 2 deletions PREP/mcip/src/xy2ll_lam.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,12 @@ SUBROUTINE xy2ll_lam (xx, yy, phi1, phi2, lambda0, phi0, phi, lambda)
term1 = DTAN ( piover4 + phi1rad/2.0d0 )
term2 = DTAN ( piover4 + phi2rad/2.0d0 )

sinphi0 = DLOG ( DCOS(phi1rad) / DCOS(phi2rad) )
sinphi0 = sinphi0 / DLOG (term2 / term1)
if ( ABS(phi1 - phi2) > 0.1) then
sinphi0 = DLOG ( DCOS(phi1rad) / DCOS(phi2rad) )
sinphi0 = sinphi0 / DLOG (term2 / term1)
else
sinphi0 = DSIN(phi1rad)
endif

sinphi0inv = 1.0d0 / sinphi0

Expand Down