jord: Geographical Position Calculations
Versions [faq]  0.1.0.0, 0.2.0.0, 0.3.0.0, 0.3.1.0, 0.4.0.0, 0.4.1.0, 0.4.2.0, 0.5.0.0, 0.6.0.0, 1.0.0.0, 2.0.0.0 

Change log  ChangeLog.md 
Dependencies  base (>=4.9 && <5), criterion, jord [details] 
License  BSD3Clause 
Copyright  2020 Cedric Liegeois 
Author  Cedric Liegeois 
Maintainer  Cedric Liegeois <ofmooseandmen@yahoo.com> 
Category  Geography 
Home page  https://github.com/ofmooseandmen/jord 
Bug tracker  https://github.com/ofmooseandmen/jord/issues 
Source repo  head: git clone https://github.com/ofmooseandmen/jord 
Uploaded  by CedricLiegeois at 20200921T07:41:34Z 
Distributions  NixOS:2.0.0.0 
Executables  jordgen, jordbenchmarks 
Downloads  3837 total (11 in the last 30 days) 
Rating  (no votes yet) [estimated by Bayesian average] 
Your Rating 

Status  Docs uploaded by user Build status unknown [no reports yet] 
Modules
[Index] [Quick Jump]
 Data
 Geo
 Jord
 Data.Geo.Jord.Angle
 Data.Geo.Jord.Duration
 Data.Geo.Jord.Ellipsoid
 Data.Geo.Jord.Ellipsoids
 Data.Geo.Jord.Geocentric
 Data.Geo.Jord.Geodesic
 Data.Geo.Jord.Geodetic
 Data.Geo.Jord.GreatCircle
 Data.Geo.Jord.Kinematics
 Data.Geo.Jord.LatLong
 Data.Geo.Jord.Length
 Data.Geo.Jord.Local
 Data.Geo.Jord.Math3d
 Data.Geo.Jord.Model
 Data.Geo.Jord.Models
 Data.Geo.Jord.Polygon
 Data.Geo.Jord.Positions
 Data.Geo.Jord.Rotation
 Data.Geo.Jord.Speed
 Data.Geo.Jord.Triangle
 Data.Geo.Jord.Tx
 Data.Geo.Jord.Txs
 Jord
 Geo
Downloads
 jord2.0.0.0.tar.gz [browse] (Cabal source package)
 Package description (as included in the package)
Readme for jord2.0.0.0
[back to package description]Jord  Geographical Position Calculations
Jord [Swedish] is Earth [English]
What is this?
Jord is a Haskell library that implements various geographical position calculations using the algorithms described in Gade, K. (2010)  A Nonsingular Horizontal Position Representation, Shudde, Rex H. (1986)  Some tactical algorithms for spherical geometry and Vincenty, T. (1975)  Direct and Inverse Solutions of Geodesics on the Ellipsoid:
 conversion between ECEF (earthcentred, earthfixed), latitude/longitude and nvector positions for spherical and ellipsoidal earth model,
 conversion between latitude/longitude and nvector positions,
 local frames (body; local level, wander azimuth; north, east, down): delta between positions, target position from reference position and delta,
 great circles: surface distance, initial & final bearing, interpolated position, great circle intersections, cross track distance, ...,
 geodesic: surface distance, initial & final bearing and destination,
 kinematics: position from p0, bearing and speed, closest point of approach between tracks, intercept (time, speed, minimum speed),
 transformation between coordinate systems (both fixed and timedependent).
 polygon triangulation
How do I build it?
If you have Stack, then:
$ stack build test
If you have Cabal, then:
$ cabal v2build
$ cabal v2test
How do I use it?
Import the core functional modules as qualified:
import qualified Data.Geo.Jord.Angle as Angle
import qualified Data.Geo.Jord.Duration as Duration
import qualified Data.Geo.Jord.Geocentric as Geocentric
import qualified Data.Geo.Jord.Geodetic as Geodetic
import qualified Data.Geo.Jord.Geodesic as Geodesic
import qualified Data.Geo.Jord.GreatCircle as GreatCircle
import qualified Data.Geo.Jord.Kinematics as Kinematics
import qualified Data.Geo.Jord.Length as Length
import qualified Data.Geo.Jord.Local as Local
import qualified Data.Geo.Jord.Polygon as Polygon
import qualified Data.Geo.Jord.Positions as Positions
import qualified Data.Geo.Jord.Speed as Speed
import qualified Data.Geo.Jord.Tx as Tx
Note: modules can be selectively imported as nonqualified, but i.m.o. the code reads better when modules are imported qualified.
Import models and transformation parameters:
import Data.Geo.Jord.Model (Epoch(..))
import Data.Geo.Jord.Models
import qualified Data.Geo.Jord.Txs as Txs
Solutions to the 10 examples from NavLab
Example 1: A and B to delta
Given two positions, A and B as latitudes, longitudes and depths relative to Earth, E.
Find the exact vector between the two positions, given in meters north, east, and down, and find the direction (azimuth) to B, relative to north. Assume WGS84 ellipsoid. The given depths are from the ellipsoid surface. Use position A to define north, east, and down directions. (Due to the curvature of Earth and different directions to the North Pole, the north, east, and down directions will change (relative to Earth) for different places. A must be outside the poles for the north and east directions to be defined.)
example1 :: IO()
example1 = do
let pA = Geodetic.latLongHeightPos 1 2 (Length.metres 3) WGS84
let pB = Geodetic.latLongHeightPos 4 5 (Length.metres 6) WGS84
let ned = Local.nedBetween pA pB
 Ned {north = 331.730863099km, east = 332.998501491km, down = 17.39830421km}
let slantRange = Local.slantRange ned
 470.357383823km
let bearing = Local.bearing ned
 45°6'33.346"
let elevation = Local.elevation ned
 2°7'11.381"
putStrLn ("NavLab, Example1: A and B to delta\n\
\ delta = " ++ (show ned) ++ "\n\
\ slantRange = " ++ (show slantRange) ++ "\n\
\ bearing = " ++ (show bearing) ++ "\n\
\ elevation = " ++ (show elevation) ++ "\n")
Example 2: B and delta to C
A radar or sonar attached to a vehicle B (Body coordinate frame) measures the distance and direction to an object C. We assume that the distance and two angles (typically bearing and elevation relative to B) are already combined to the vector p_BC_B (i.e. the vector from B to C, decomposed in B). The position of B is given as n_EB_E and z_EB, and the orientation (attitude) of B is given as R_NB (this rotation matrix can be found from roll/pitch/yaw by using zyx2R).
Find the exact position of object C as nvector and depth ( n_EC_E and z_EC ), assuming Earth ellipsoid with semimajor axis a and flattening f. For WGS72, use a = 6 378 135 m and f = 1/298.26.
example2 :: IO()
example2 = do
let frameB =
Local.frameB
(Angle.decimalDegrees 40)
(Angle.decimalDegrees 20)
(Angle.decimalDegrees 30)
let pB = Geodetic.nvectorHeightPos 1 2 3 (Length.metres 400) WGS72
let delta = Local.deltaMetres 3000 2000 100
let pC = Local.destination pB frameB delta
 53°18'46.839"N,63°29'6.179"E 406.006017m (WGS72)
putStrLn ("NavLab, Example 2: B and delta to C\n\
\ pC = " ++ (show pC) ++ "\n")
Example 3: ECEFvector to geodetic latitude
Position B is given as an “ECEFvector” p_EB_E (i.e. a vector from E, the center of the Earth, to B, decomposed in E). Find the geodetic latitude, longitude and height (latEB, lonEB and hEB), assuming WGS84 ellipsoid.
example3 :: IO()
example3 = do
let ecef = Geocentric.metresPos 5733900.0 (6371000.0) 7008100.000000001 WGS84
let geod = Positions.toGeodetic ecef
 39°22'43.495"N,48°0'46.035"W 4702.059834295km (WGS84)
putStrLn ("NavLab, 3: ECEFvector to geodetic latitude\n\
\ geodetic pos = " ++ (show geod) ++ "\n")
Example 4: Geodetic latitude to ECEFvector
Geodetic latitude, longitude and height are given for position B as latEB, lonEB and hEB, find the ECEFvector for this position, p_EB_E.
example4 :: IO()
example4 = do
let geod = Geodetic.latLongHeightPos 1 2 (Length.metres 3) WGS84
let ecef = Positions.toGeocentric geod
 Position {gx = 6373.290277218km, gy = 222.560200675km, gz = 110.568827182km, model = WGS84}
putStrLn ("NavLab, 4: Geodetic latitude to ECEFvector\n\
\ geocentric pos = " ++ (show ecef) ++ "\n")
Example 5: Surface distance
Find the surface distance sAB (i.e. great circle distance) between two positions A and B. The heights of A and B are ignored, i.e. if they don’t have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B. The Euclidean distance (chord length) dAB should also be found. Use Earth radius 6371e3 m. Compare the results with exact calculations for the WGS84 ellipsoid.
example5 :: IO()
example5 = do
let pA = Geodetic.s84Pos 88 0
let pB = Geodetic.s84Pos 89 (170)
let distance = GreatCircle.distance pA pB
 332.456901835km
putStrLn ("NavLab, 5: Surface distance\n\
\ distance = " ++ (show distance) ++ "\n")
Example 6: Interpolated position
Given the position of B at time t0 and t1, n_EB_E(t0) and n_EB_E(t1).
Find an interpolated position at time ti, n_EB_E(ti). All positions are given as nvectors.
example6 :: IO()
example6 = do
let pA = Geodetic.s84Pos 89 0
let pB = Geodetic.s84Pos 89 180
let f = 0.6
let interpolated = GreatCircle.interpolated pA pB f
 89°47'59.929"N,180°0'0.000"E (S84)
putStrLn ("NavLab, 6: Interpolated position\n\
\ interpolated = " ++ (show interpolated) ++ "\n")
Example 7: Mean position
Three positions A, B, and C are given as nvectors n_EA_E, n_EB_E, and n_EC_E. Find the mean position, M, given as n_EM_E. Note that the calculation is independent of the depths of the positions.
example7 :: IO()
example7 = do
let ps =
[ Geodetic.s84Pos 90 0
, Geodetic.s84Pos 60 10
, Geodetic.s84Pos 50 (20)
]
let mean = GreatCircle.mean ps
 Just 67°14'10.150"N,6°55'3.040"W (S84)
putStrLn ("NavLab, Example 7: Mean position\n\
\ mean = " ++ (show mean) ++ "\n")
Example 8: A and azimuth/distance to B
We have an initial position A, direction of travel given as an azimuth (bearing) relative to north (clockwise), and finally the distance to travel along a great circle given as sAB. Use Earth radius 6371e3 m to find the destination point B.
In geodesy this is known as “The first geodetic problem” or “The direct geodetic problem” for a sphere, and we see that this is similar to Example 2, but now the delta is given as an azimuth and a great circle distance. (“The second/inverse geodetic problem” for a sphere is already solved in Examples 1 and 5.)
example8 :: IO()
example8 = do
let p = Geodetic.s84Pos 80 (90)
let bearing = Angle.decimalDegrees 200
let distance = Length.metres 1000
let dest = GreatCircle.destination p bearing distance
 79°59'29.575"N,90°1'3.714"W (S84)
putStrLn ("NavLab, Example 8: A and azimuth/distance to B\n\
\ destination = " ++ (show dest) ++ "\n")
Example 9: Intersection of two paths
Define a path from two given positions (at the surface of a spherical Earth), as the great circle that goes through the two points.
Path A is given by A1 and A2, while path B is given by B1 and B2.
Find the position C where the two great circles intersect.
example9 :: IO()
example9 = do
let a1 = Geodetic.s84Pos 51.885 0.235
let a2 = Geodetic.s84Pos 48.269 13.093
let b1 = Geodetic.s84Pos 49.008 2.549
let b2 = Geodetic.s84Pos 56.283 11.304
let ga = GreatCircle.through a1 a2
let gb = GreatCircle.through b1 b2
let intersections = GreatCircle.intersections <$> ga <*> gb
 Just (Just (50°54'6.260"N,4°29'39.052"E (S84),50°54'6.260"S,175°30'20.947"W (S84)))
let ma = GreatCircle.minorArc a1 a2
let mb = GreatCircle.minorArc b1 b2
let intersection = GreatCircle.intersection <$> ma <*> mb
 Just (Just 50°54'6.260"N,4°29'39.052"E (S84))
putStrLn ("NavLab, Example 9: Intersection of two paths\n\
\ great circle intersections = " ++ (show intersections) ++ "\n\
\ minor arc intersection = " ++ (show intersection) ++ "\n")
Example 10: Cross track distance
*Path A is given by the two positions A1 and A2 (similar to the previous example).
Find the cross track distance sxt between the path A (i.e. the great circle through A1 and A2) and the position B (i.e. the shortest distance at the surface, between the great circle and B).
example10 :: IO()
example10 = do
let p = Geodetic.s84Pos 1 0.1
let gc = GreatCircle.through (Geodetic.s84Pos 0 0) (Geodetic.s84Pos 10 0)
let sxt = fmap (\g > GreatCircle.crossTrackDistance p g) gc
 Just 11.117814411km
putStrLn ("NavLab, Example 10: Cross track distance\n\
\ cross track distance = " ++ (show sxt) ++ "\n")
Solutions to the geodesic problems (Vincenty)
Example 11: Inverse problem
Given the coordinates of the two points (Φ1, L1) and (Φ2, L2), the inverse problem finds the azimuths α1, α2 and the ellipsoidal distance s.
example11 :: IO()
example11 = do
let pA = Geodetic.wgs84Pos 88 0
let pB = Geodetic.wgs84Pos 89 (170)
let inv = Geodesic.inverse pA pB
let initialBearing = fmap Geodesic.initialBearing inv
 Just (Just 356°40'8.701")
let finalBearing = fmap Geodesic.finalBearing inv
 Just (Just 186°40'19.615")
let distance = fmap Geodesic.length inv
 Just 333.947509469km
putStrLn ("Geodesic, Example 11: Inverse problem\n\
\ initial bearing = " ++ (show initialBearing) ++ "\n\
\ final bearing = " ++ (show finalBearing) ++ "\n\
\ distance = " ++ (show distance) ++ "\n")
Example 12: Direct problem
Given an initial point (Φ1, L1) and initial azimuth, α1, and a distance, s, along the geodesic the problem is to find the end point (Φ2, L2) and azimuth, α2.
example12 :: IO()
example12 = do
let p = Geodetic.wgs84Pos 80 (90)
let bearing = Angle.decimalDegrees 200
let distance = Length.metres 1000
let dct = Geodesic.direct p bearing distance
let destination = fmap Geodesic.endPosition dct
 Just 79°59'29.701"N,90°1'3.436"W (WGS84)
let finalBearing = fmap Geodesic.finalBearing dct
 Just (Just 199°58'57.528")
putStrLn ("Geodesic, Example 12: Direct problem\n\
\ destination = " ++ (show destination) ++ "\n\
\ final bearing = " ++ (show finalBearing) ++ "\n")
Solutions to kinematics problems
Example 13: Closest point of approach
The Closest Point of Approach (CPA) refers to the positions at which two dynamically moving objects reach their closest possible distance.
example13 :: IO()
example13 = do
let ownship = Kinematics.Track
(Geodetic.s84Pos 20 (60))
(Angle.decimalDegrees 10)
(Speed.knots 15)
let intruder = Kinematics.Track
(Geodetic.s84Pos 34 (50))
(Angle.decimalDegrees 220)
(Speed.knots 300)
let cpa = Kinematics.cpa ownship intruder
let timeToCpa = fmap Kinematics.timeToCpa cpa
 Just 3H9M56.155S
let distanceAtCpa = fmap Kinematics.distanceAtCpa cpa
 Just 124.231730834km
let cpaOwnshipPosition = fmap Kinematics.cpaOwnshipPosition cpa
 Just 20°46'43.641"N,59°51'11.225"W (S84)
let cpaIntruderPosition = fmap Kinematics.cpaIntruderPosition cpa
 Just 21°24'8.523"N,60°50'48.159"W (S84)
putStrLn ("Kinematics, Example 13: Closest point of approach\n\
\ time to CPA = " ++ (show timeToCpa) ++ "\n\
\ distance at CPA = " ++ (show distanceAtCpa) ++ "\n\
\ CPA ownship pos = " ++ (show cpaOwnshipPosition) ++ "\n\
\ CPA intruder pos = " ++ (show cpaIntruderPosition) ++ "\n")
Example 14: Speed required to intercept target
Inputs are the initial latitude and longitude of an interceptor and a target, and the target course and speed. Also input is the time of the desired intercept. Outputs are the speed required of the interceptor, the course of the interceptor, the distance travelled to intercept, and the latitude and longitude of the intercept.
example14 :: IO()
example14 = do
let track = Kinematics.Track
(Geodetic.s84Pos 34 (50))
(Angle.decimalDegrees 220)
(Speed.knots 600)
let interceptor = Geodetic.s84Pos 20 (60)
let interceptTime = Duration.seconds 2700
let intercept = Kinematics.interceptByTime track interceptor interceptTime
let distanceToIntercept = fmap Kinematics.distanceToIntercept intercept
 Just 1015.302358852km
let interceptPosition = fmap Kinematics.interceptPosition intercept
 Just 28°8'12.046"N,55°27'21.411"W (S84)
let interceptorBearing = fmap Kinematics.interceptorBearing intercept
 Just 26°7'11.649"
let interceptorSpeed = fmap Kinematics.interceptorSpeed intercept
 Just 1353.736478km/h
putStrLn ("Kinematics, Example 14: Speed required to intercept target\n\
\ distance to intercept = " ++ (show distanceToIntercept) ++ "\n\
\ intercept position = " ++ (show interceptPosition) ++ "\n\
\ interceptor bearing = " ++ (show interceptorBearing) ++ "\n\
\ interceptor speed = " ++ (show interceptorSpeed) ++ "\n")
Example 15: Time required to intercept target
Inputs are the initial latitude and longitude of an interceptor and a target, and the target course and speed. For a given interceptor speed, it may or may not be possible to make an intercept.
The first algorithm is to compute the minimum interceptor speed required to achieve intercept and the time required to make such and intercept.
The second algorithm queries the user to input an interceptor speed. If the speed is at least that required for intercept then the time required to intercept is computed.
example15 :: IO()
example15 = do
let track = Kinematics.Track
(Geodetic.s84Pos 34 (50))
(Angle.decimalDegrees 220)
(Speed.knots 600)
let interceptor = Geodetic.s84Pos 20 (60)
let minIntercept = Kinematics.intercept track interceptor
let minTimeToIntercept = fmap Kinematics.timeToIntercept minIntercept
 Just 1H39M53.831S
let minDistanceToIntercept = fmap Kinematics.distanceToIntercept minIntercept
 Just 162.294627463km
let minInterceptPosition = fmap Kinematics.interceptPosition minIntercept
 Just 20°43'42.305"N,61°20'56.848"W (S84)
let minInterceptorBearing = fmap Kinematics.interceptorBearing minIntercept
 Just 300°10'18.053"
let minInterceptorSpeed = fmap Kinematics.interceptorSpeed minIntercept
 Just 97.476999km/h
putStrLn ("Kinematics, Example 15: Time required to intercept target (min)\n\
\ time to intercept = " ++ (show minTimeToIntercept) ++ "\n\
\ distance to intercept = " ++ (show minDistanceToIntercept) ++ "\n\
\ intercept position = " ++ (show minInterceptPosition) ++ "\n\
\ interceptor bearing = " ++ (show minInterceptorBearing) ++ "\n\
\ interceptor speed = " ++ (show minInterceptorSpeed) ++ "\n")
let interceptSpeed = Speed.knots 700
let intercept = Kinematics.interceptBySpeed track interceptor interceptSpeed
let timeToIntercept = fmap Kinematics.timeToIntercept intercept
 Just 0H46M4.692S
let distanceToIntercept = fmap Kinematics.distanceToIntercept intercept
 Just 995.596069189km
let interceptPosition = fmap Kinematics.interceptPosition intercept
 Just 27°59'36.764"N,55°34'43.852"W(S84)
let interceptorBearing = fmap Kinematics.interceptorBearing intercept
 Just 25°56'7.484"
putStrLn ("Kinematics, Example 15: Time required to intercept target\n\
\ time to intercept = " ++ (show timeToIntercept) ++ "\n\
\ distance to intercept = " ++ (show distanceToIntercept) ++ "\n\
\ intercept position = " ++ (show interceptPosition) ++ "\n\
\ interceptor bearing = " ++ (show interceptorBearing) ++ "\n")
Solutions to coordinates transformation problems
Example 16: Transformation between fixed models
Convert the coordinates of a geocentric position between the WGS84 model and the NAD83 model.
example16 :: IO()
example16 = do
let pWGS84 = Geocentric.metresPos 4193790.895437 454436.195118 4768166.813801 WGS84
 using explicit parameters:
let tx = Tx.params Txs.from_WGS84_to_NAD83
let pNAD83 = Positions.transform' pWGS84 NAD83 tx
 Position {gx = 4193.792080781km, gy = 454.433921298km, gz = 4768.16615479km, model = NAD83}
 using the transformation graph:
let pNAD83' = Positions.transform pWGS84 NAD83 Txs.fixed
 Just (Position {gx = 4193.792080781km, gy = 454.433921298km, gz = 4768.16615479km, model = NAD83})
putStrLn ("Coordinates transformation, Example 16: WGS84 > NAD83\n\
\ WGS84 = " ++ (show pWGS84) ++ "\n\
\ NAD83 = " ++ (show pNAD83) ++ "\n\
\ NAD83 = " ++ (show pNAD83') ++ "\n")
Example 17: Transformation between time dependent models
Convert the coordinates of a geocentric position between the ITRF2014, ITRF2000 and NAD83 (CORS96) models
example17 :: IO()
example17 = do
let pITRF2014 = Geocentric.metresPos 4193790.895437 454436.195118 4768166.813801 ITRF2014
 using explicit parameters:
let tx = Tx.params Txs.from_ITRF2014_to_ITRF2000
let pITRF2000 = Positions.transformAt' pITRF2014 (Epoch 2014.0) ITRF2000 tx
 Position {gx = 4193.790907273km, gy = 454.436197881km, gz = 4768.166792308km, model = ITRF2000}
 using the transformation graph (goes via ITRF2000):
let pNAD83 = Positions.transformAt pITRF2014 (Epoch 2014.0) NAD83_CORS96 Txs.timeDependent
 Just (Position {gx = 4193.791801305km, gy = 454.433876376km, gz = 4768.166397281km, model = NAD83_CORS96})
putStrLn ("Coordinates transformation, Example 17: ITRF2014 > ITRF2000 > NAD83 (CORS96)\n\
\ ITRF2014 = " ++ (show pITRF2014) ++ "\n\
\ ITRF2000 = " ++ (show pITRF2000) ++ "\n\
\ NAD83 (CORS96) = " ++ (show pNAD83) ++ "\n")
Solutions to polygon triangulation
Example 18: Triangulation of a simple polygon
example18 :: IO()
example18 = do
let p = Polygon.simple
[ Geodetic.s84Pos 45 45
, Geodetic.s84Pos (45) 45
, Geodetic.s84Pos (45) (45)
, Geodetic.s84Pos 45 (45)
]
let ts = fmap Polygon.triangulate p
 Right [ Triangle 45°0'0.000"N,45°0'0.000"W (S84) 45°0'0.000"N,45°0'0.000"E (S84) 45°0'0.000"S,45°0'0.000"E (S84)
 , Triangle 45°0'0.000"S,45°0'0.000"E (S84) 45°0'0.000"S,45°0'0.000"W (S84) 45°0'0.000"N,45°0'0.000"W (S84)
 ]
putStrLn ("Polygon triangulation, Example 18: Simple polygon\n\
\ triangles = " ++ (show ts) ++ "\n")
Example 19: Triangulation of a circle
example19 :: IO()
example19 = do
let p = Polygon.circle (Geodetic.s84Pos 0 0) (Length.kilometres 10) 10
let ts = fmap Polygon.triangulate p
 Right [ Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°5'23.755"N,0°0'0.000"E (S84) 0°4'21.923"N,0°3'10.298"E (S84)
 , Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°4'21.923"N,0°3'10.298"E (S84) 0°1'40.045"N,0°5'7.909"E (S84)
 , Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°1'40.045"N,0°5'7.909"E (S84) 0°1'40.045"S,0°5'7.909"E (S84)
 , Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°1'40.045"S,0°5'7.909"E (S84) 0°4'21.923"S,0°3'10.298"E (S84)
 , Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°4'21.923"S,0°3'10.298"E (S84) 0°5'23.755"S,0°0'0.000"E (S84)
 , Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°5'23.755"S,0°0'0.000"E (S84) 0°4'21.923"S,0°3'10.298"W (S84)
 , Triangle 0°4'21.923"N,0°3'10.298"W (S84) 0°4'21.923"S,0°3'10.298"W (S84) 0°1'40.045"S,0°5'7.909"W (S84)
 , Triangle 0°1'40.045"S,0°5'7.909"W (S84) 0°1'40.045"N,0°5'7.909"W (S84) 0°4'21.923"N,0°3'10.298"W (S84)
 ]
putStrLn ("Polygon triangulation, Example 19: Circle\n\
\ triangles = " ++ (show ts) ++ "\n")
Example 20: Triangulation of an arc
example20 :: IO()
example20 = do
let p = Polygon.arc (Geodetic.s84Pos 0 0) (Length.kilometres 10) (Angle.decimalDegrees 10) (Angle.decimalDegrees 20) 5
let ts = fmap Polygon.triangulate p
 Right [ Triangle 0°5'4.230"N,0°1'50.730"E (S84) 0°5'18.836"N,0°0'56.219"E (S84) 0°5'16.081"N,0°1'10.073"E (S84)
 , Triangle 0°5'4.230"N,0°1'50.730"E (S84) 0°5'16.081"N,0°1'10.073"E (S84) 0°5'12.723"N,0°1'23.794"E (S84)
 , Triangle 0°5'12.723"N,0°1'23.794"E (S84) 0°5'8.770"N,0°1'37.355"E (S84) 0°5'4.230"N,0°1'50.730"E (S84)
 ]
putStrLn ("Polygon triangulation, Example 20: Arc\n\
\ triangles = " ++ (show ts) ++ "\n")
main :: IO()
main = do
example1
example2
example3
example4
example5
example6
example7
example8
example9
example10
example11
example12
example13
example14
example15
example16
example17
example18
example19
example20