shifts geodetic coordinates relative to a given source datum to geodetic coordinates relative to WGS84.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(Coordinate), | intent(in) | :: | input | |||
real(kind=float), | intent(in) | :: | latIn | |||
real(kind=float), | intent(in) | :: | lonIn | |||
real(kind=float), | intent(out) | :: | WGS84lat | |||
real(kind=float), | intent(out) | :: | WGS84lon |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
real(kind=float), | public, | parameter | :: | MOLODENSKY_MAX | = | 89.75*degToRad |
Polar limit |
real(kind=float), | public | :: | WGS84_a |
Semi-major axis of WGS84 ellipsoid in meters |
|||
real(kind=float), | public | :: | WGS84_f |
Flattening of WGS84 ellisoid |
|||
real(kind=float), | public | :: | WGS84hgt | ||||
real(kind=float), | public | :: | a |
Semi-major axis of ellipsoid in meters |
|||
real(kind=float), | public | :: | da |
Difference in semi-major axes |
|||
real(kind=float), | public | :: | df |
Difference in flattening |
|||
real(kind=float), | public | :: | dx | ||||
real(kind=float), | public | :: | dy | ||||
real(kind=float), | public | :: | dz | ||||
real(kind=float), | public | :: | f |
Flattening of ellipsoid |
|||
real(kind=float), | public | :: | hgtIn |
SUBROUTINE GeodeticShiftToWGS84 & ! (input, latIn, lonIn, WGS84lat, WGS84lon) IMPLICIT NONE ! Arguments with intent(in): TYPE (Coordinate), INTENT(IN) :: input REAL (KIND = float), INTENT(IN) :: latIn REAL (KIND = float), INTENT(IN) :: lonIn ! Arguments with intent(out): REAL (KIND = float), INTENT(OUT) :: WGS84lat, WGS84lon ! Local variables: REAL (KIND = float) :: hgtIn REAL (KIND = float) :: WGS84_a !! Semi-major axis of WGS84 ellipsoid in meters REAL (KIND = float) :: WGS84_f !! Flattening of WGS84 ellisoid REAL (KIND = float) :: a !!Semi-major axis of ellipsoid in meters REAL (KIND = float) :: f !!Flattening of ellipsoid REAL (KIND = float) :: da !!Difference in semi-major axes REAL (KIND = float) :: df !!Difference in flattening REAL (KIND = float) :: dx REAL (KIND = float) :: dy REAL (KIND = float) :: dz REAL (KIND = float) :: WGS84hgt ! Local parameters: REAL (KIND = float), PARAMETER :: MOLODENSKY_MAX = 89.75 * degToRad !! Polar limit !------------end of declaration------------------------------------------------ !initialize hgtIn = 0. WGS84_a = 6378137.0 WGS84_f = 1. / 298.257223563 a = input % system % ellipsoid % a f = input % system % ellipsoid % f !If input datum is WGS84, just copy to output IF ( input % system % datum % code == WGS84 ) THEN WGS84lat = latIn WGS84lon = lonIn RETURN END IF !If latitude is within limits, apply Molodensky IF ( latIn <= MOLODENSKY_MAX .AND. & latIn >= - MOLODENSKY_MAX ) THEN da = WGS84_a - a df = WGS84_f - f dx = input % system % datum % dx dy = input % system % datum % dy dz = input % system % datum % dz CALL MolodenskyShift(WGS84_a, da, WGS84_f, df, dx, dy, dz, latIn, lonIn, & hgtIn, WGS84lat, WGS84lon, WGS84hgt) ELSE !apply 3 steps method WRITE (*,*) 'Latitude is outside Molodensky limits' WRITE (*,*) '3-step datum transformation not yet implemented' STOP END IF END SUBROUTINE GeodeticShiftToWGS84