2 people like it.

Azimuthal equidistant projection with measures

Simple version of the azimuthal equidistant projection (see also http://fssnip.net/lA) but with measures. This avoids mixing of degrees and radians and longitudes/x and latitudes/y

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
    open System

    [<Measure>] type deg
    [<Measure>] type rad
    [<Measure>] type x
    [<Measure>] type y
    type lon = float<x deg>
    type lat = float<y deg>
    let inline degToRad (d:float<'u deg>) = d*0.0174532925199433<rad/deg>
    let cos (d:float<'u rad>) = Math.Cos(float d)
    let sin (d:float<'u rad>) = Math.Sin(float d)
    
    let project (centerlon:lon) (centerlat:lat) (lon:lon) (lat:lat) =
        // http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
        // http://www.radicalcartography.net/?projectionref
        let t = degToRad lat
        let l = degToRad lon
        let t1 = degToRad centerlat // latitude center of projection
        let l0 = degToRad centerlon // longitude center of projection
        let c = Math.Acos ((sin t1) * (sin t) + (cos t1) * (cos t) * (cos (l-l0)))
        let k = c / (sin c)
        let x:float<x> = k * (cos t) * (sin (l-l0)) 
                         |> LanguagePrimitives.FloatWithMeasure
        let y:float<y> = k * (cos t1) * (sin t) - (sin t1) * (cos t) * (cos (l-l0)) 
                         |> LanguagePrimitives.FloatWithMeasure
        (x, y)
namespace System
Multiple items
type MeasureAttribute =
  inherit Attribute
  new : unit -> MeasureAttribute

Full name: Microsoft.FSharp.Core.MeasureAttribute

--------------------
new : unit -> MeasureAttribute
[<Measure>]
type deg

Full name: Script.deg
[<Measure>]
type rad

Full name: Script.rad
[<Measure>]
type x

Full name: Script.x
[<Measure>]
type y

Full name: Script.y
type lon = float<deg x>

Full name: Script.lon
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
type lat = float<deg y>

Full name: Script.lat
val degToRad : d:float<'u> -> float<'u rad/deg>

Full name: Script.degToRad
val d : float<'u>
val cos : d:float<'u> -> float

Full name: Script.cos
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
Math.Cos(d: float) : float
val sin : d:float<'u> -> float

Full name: Script.sin
Math.Sin(a: float) : float
val project : centerlon:lon -> centerlat:lat -> lon:lon -> lat:lat -> float<x> * float<y>

Full name: Script.project
val centerlon : lon
val centerlat : lat
Multiple items
val lon : lon

--------------------
type lon = float<deg x>

Full name: Script.lon
Multiple items
val lat : lat

--------------------
type lat = float<deg y>

Full name: Script.lat
val t : float<rad y>
val l : float<rad x>
val t1 : float<rad y>
val l0 : float<rad x>
val c : float
Math.Acos(d: float) : float
val k : float
val x : float<x>
module LanguagePrimitives

from Microsoft.FSharp.Core
val FloatWithMeasure : float -> float<'Measure>

Full name: Microsoft.FSharp.Core.LanguagePrimitives.FloatWithMeasure
val y : float<y>
Raw view Test code New version

More information

Link:http://fssnip.net/lD
Posted:10 years ago
Author:Samuel Bosch
Tags: math , spatial