ConvertGeodeticToSwiss Subroutine

public subroutine ConvertGeodeticToSwiss(lon, lat, k, lon0, lat0, azimuth, a, e, falseN, falseE, x, y)

Uses

The subroutine converts geodetic (latitude and longitude) to Swiss Oblique Cylindrical projection(easting and northing) coordinates

Reference:

Federal Office of Topography, Formulas and constants for the calculation of the Swiss conformal cylindrical projection and for the transormation between coordinate systems http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys.html

Arguments

Type IntentOptional Attributes Name
real(kind=float), intent(in) :: lon

geodetic longitude [radians]
real(kind=float), intent(in) :: lat

geodetic latitude [radians]
real(kind=float), intent(in) :: k

scale factor
real(kind=float), intent(in) :: lon0

longitude of center [radians]
real(kind=float), intent(in) :: lat0

latitude of center [radians]
real(kind=float), intent(in) :: azimuth

azimuth of centerline
real(kind=float), intent(in) :: a

semimajor axis [m]
real(kind=float), intent(in) :: e

eccentricity
real(kind=float), intent(in) :: falseN

false northing
real(kind=float), intent(in) :: falseE

false easting
real(kind=float), intent(out) :: x

easting coordinate [m]
real(kind=float), intent(out) :: y

northing coordinate [m]

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: K0
real(kind=float), public :: R
real(kind=float), public :: S
real(kind=float), public :: alpha
real(kind=float), public :: b
real(kind=float), public :: b0
real(kind=float), public :: b1
real(kind=float), public :: e2
integer(kind=short), public :: i
real(kind=float), public :: l
real(kind=float), public :: l1
real(kind=float), public :: xbig
real(kind=float), public :: ybig

Source Code

SUBROUTINE ConvertGeodeticToSwiss &
!!
(lon, lat, k, lon0, lat0, azimuth, a, e, falseN, falseE, x, y)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString


IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (IN) :: lat !!geodetic latitude [radians]
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: lon0 !!longitude of center [radians]
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of center [radians]
REAL (KIND = float), INTENT (IN) :: azimuth !!azimuth of centerline
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: falseN !!false northing
REAL (KIND = float), INTENT (IN) :: falseE !!false easting


!Arguments with intent (out):
REAL (KIND = float), INTENT (OUT) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (OUT) :: y !!northing coordinate [m]

!Local declarations:
REAL (KIND = float) :: e2
REAL (KIND = float) :: R
REAL (KIND = float) :: alpha
REAL (KIND = float) :: b0
REAL (KIND = float) :: K0
REAL (KIND = float) :: xbig, ybig
REAL (KIND = float) :: l1
REAL (KIND = float) :: b1
REAL (KIND = float) :: l
REAL (KIND = float) :: b
REAL (KIND = float) :: S
INTEGER (KIND = short) :: i

!------------end of declaration------------------------------------------------

!eccentricity squared
e2 = e * e

!Radius of the projection sphere
R = a * (1. - e2)**0.5 / (1. - e2 * (SIN(lat0))**2.)

!Relat. between longitude on sphere and on ellipsoid
alpha = ( 1. + e2 / (1. - e2) * COS(lat0)**4. )**0.5

!Latitude of the fundamental point on the sphere
b0 = ASIN( SIN(lat0)/alpha )

!Constant of the latitude
K0 = LOG(TAN(pi/4. + b0/2.)) - alpha * LOG(TAN(pi/4. + lat0/2.)) + &
     alpha * e / 2. * LOG( (1. + e * SIN(lat0)) / (1. - e * SIN(lat0)) )


!ellipsoid to sphere
!auxiliary value
S = alpha * LOG(TAN(pi/4. + lat/2.)) - alpha * e / 2. * &
    LOG( (1. + e * SIN(lat)) / (1. - e * SIN(lat)) ) + K0
!spherical latitude
b = 2. * ( ATAN(EXP(S)) - pi/4. )
!spherical longitude
l = alpha * (lon - lon0)

!equator system -> pseudo-equator system
l1 = ATAN( SIN(l) / (SIN(b0) * TAN(b) + COS(b0) * COS(l)) )

b1 = ASIN(COS(b0) * SIN(b) - SIn(b0) * COS(b) * COS(l))

!sphere -> projection plane
xbig = R * l1
ybig = R / 2. * LOG((1. + SIN(b1)) / (1. - SIN(b1)))

x = xbig + falseE
y = ybig + falseN


END SUBROUTINE ConvertGeodeticToSwiss