diff --git a/changelog b/changelog
new file mode 100644
--- /dev/null
+++ b/changelog
@@ -0,0 +1,9 @@
+0.1.3
+
+  * Update lens dependency to >= 4.0
+
+  * Update coordinate dependency to >= 0.0.16
+  
+0.1.2
+
+* Release under NICTA open source.
diff --git a/geodetic.cabal b/geodetic.cabal
--- a/geodetic.cabal
+++ b/geodetic.cabal
@@ -1,5 +1,5 @@
 name:               geodetic
-version:            0.1.2
+version:            0.1.3
 license:            BSD3
 license-File:       etc/LICENCE
 author:             Tony Morris <ʇǝu˙sıɹɹoɯʇ@ןןǝʞsɐɥ>
@@ -15,6 +15,7 @@
 bug-reports:        https://github.com/NICTA/geodetic/issues
 cabal-version:      >= 1.10
 build-type:         Custom
+extra-source-files: changelog
 
 source-repository   head
   type:             git
@@ -29,8 +30,8 @@
 
   build-depends:
                     base < 5 && >= 3
-                    , lens >= 3.10
-                    , coordinate >= 0.0.8
+                    , lens >= 4.0
+                    , coordinate >= 0.0.16
                     , optional >= 0.0.1
 
   ghc-options:
@@ -68,7 +69,8 @@
                     doctest >= 0.9.7,
                     filepath >= 1.3,
                     directory >= 1.1,
-                    QuickCheck >= 2.0
+                    QuickCheck >= 2.0,
+                    template-haskell >= 2.8
 
   ghc-options:
                     -Wall
diff --git a/src/Data/Geo/Geodetic/Azimuth.hs b/src/Data/Geo/Geodetic/Azimuth.hs
--- a/src/Data/Geo/Geodetic/Azimuth.hs
+++ b/src/Data/Geo/Geodetic/Azimuth.hs
@@ -1,14 +1,18 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+
+
 -- | An azimuth in degrees between 0 and 360.
 module Data.Geo.Geodetic.Azimuth(
   Azimuth
-, HasAzimuth(..)
+, AsAzimuth(..)
 , modAzimuth
-, nAzimuth
 ) where
 
+import Control.Applicative(Applicative)
 import Prelude(Double, Bool(..), Eq, Show(..), Ord(..), id, (&&), (++), showParen, showString)
 import Data.Maybe(Maybe(..))
-import Control.Lens(Prism', Lens', prism')
+import Control.Lens(Choice, Optic', prism')
 import Text.Printf(printf)
 import Data.Fixed(mod')
 
@@ -50,34 +54,33 @@
 modAzimuth x =
   Azimuth (x `mod'` 360)
 
+class AsAzimuth p f s where
+  _Azimuth ::
+    Optic' p f s Azimuth
+
+instance AsAzimuth p f Azimuth where
+  _Azimuth =
+    id
+
 -- | A prism on azimuth to an integer between 0 and 359 inclusive.
 --
--- >>> 7 ^? nAzimuth
+-- >> 7 ^? _Azimuth
 -- Just (Azimuth 7.0000)
 --
--- >>> 0 ^? nAzimuth
+-- >> 0 ^? _Azimuth
 -- Just (Azimuth 0.0000)
 --
--- >>> 359.999 ^? nAzimuth
+-- >> 359.999 ^? _Azimuth
 -- Just (Azimuth 359.9990)
 --
--- >>> 360 ^? nAzimuth
+-- >> 360 ^? _Azimuth
 -- Nothing
 --
--- prop> all (\m -> nAzimuth # m == n) (n ^? nAzimuth)
-nAzimuth ::
-  Prism' Double Azimuth
-nAzimuth =
-  prism'
-    (\(Azimuth i) -> i)
-    (\i -> case i >= 0 && i < 360 of
-             True -> Just (Azimuth i)
-             False -> Nothing)
-
-class HasAzimuth t where
-  azimuth ::
-    Lens' t Azimuth
-
-instance HasAzimuth Azimuth where
-  azimuth =
-    id
+-- prlop> all (\m -> _Azimuth # m == n) (n ^? _Azimuth)
+instance (Choice p, Applicative f) => AsAzimuth p f Double where
+  _Azimuth =
+    prism'
+      (\(Azimuth i) -> i)
+      (\i -> case i >= 0 && i < 360 of
+               True -> Just (Azimuth i)
+               False -> Nothing)        
diff --git a/src/Data/Geo/Geodetic/Bearing.hs b/src/Data/Geo/Geodetic/Bearing.hs
--- a/src/Data/Geo/Geodetic/Bearing.hs
+++ b/src/Data/Geo/Geodetic/Bearing.hs
@@ -1,15 +1,18 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+
 -- | A bearing in degrees between 0 and 360.
 module Data.Geo.Geodetic.Bearing(
   Bearing
-, HasBearing(..)
+, AsBearing(..)
 , modBearing
-, degreeBearing
 , radianBearing
 ) where
 
+import Control.Applicative(Applicative)
 import Prelude(Double, Bool(..), Eq, Show(..), Num(..), Fractional(..), Ord(..), id, (&&), (++), (.), showString, showParen, pi)
 import Data.Maybe(Maybe(..))
-import Control.Lens(Prism', Lens', prism', iso)
+import Control.Lens(Choice, Optic', Prism', prism', iso)
 import Text.Printf(printf)
 import Data.Fixed(mod')
 
@@ -54,33 +57,6 @@
 modBearing x =
   Bearing (x `mod'` 360)
 
--- | A prism on bearing to a double between 0 inclusive and 360 exclusive.
---
--- >>> 7 ^? degreeBearing
--- Just (Bearing 7.0000)
---
--- >>> 0 ^? degreeBearing
--- Just (Bearing 0.0000)
---
--- >>> 359 ^? degreeBearing
--- Just (Bearing 359.0000)
---
--- >>> 359.997 ^? degreeBearing
--- Just (Bearing 359.9970)
---
--- >>> 360 ^? degreeBearing
--- Nothing
---
--- prop> all (\m -> degreeBearing # m == n) (n ^? degreeBearing)
-degreeBearing ::
-  Prism' Double Bearing
-degreeBearing =
-  prism'
-    (\(Bearing i) -> i)
-    (\i -> case i >= 0 && i < 360 of
-             True -> Just (Bearing i)
-             False -> Nothing)
-
 -- | A prism on bearing to a double between 0 and π exclusive.
 --
 -- >>> (2 * pi - 0.0000000001) ^? radianBearing
@@ -106,12 +82,39 @@
 radianBearing ::
   Prism' Double Bearing
 radianBearing =
-  iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . degreeBearing
+  iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Bearing
 
-class HasBearing t where
-  bearing ::
-    Lens' t Bearing
+class AsBearing p f s where
+  _Bearing ::
+    Optic' p f s Bearing
 
-instance HasBearing Bearing where
-  bearing =
+instance AsBearing p f Bearing where
+  _Bearing =
     id
+
+-- | A prism on bearing to a double between 0 inclusive and 360 exclusive.
+--
+-- >>> (7 :: Double) ^? _Bearing
+-- Just (Bearing 7.0000)
+--
+-- >>> (0 :: Double) ^? _Bearing
+-- Just (Bearing 0.0000)
+--
+-- >>> (359 :: Double) ^? _Bearing
+-- Just (Bearing 359.0000)
+--
+-- >>> (359.997 :: Double) ^? _Bearing
+-- Just (Bearing 359.9970)
+--
+-- >>> (360 :: Double) ^? _Bearing
+-- Nothing
+--
+-- prop> all (\m -> _Bearing # m == n) ((n :: Double) ^? _Bearing)
+instance (Choice p, Applicative f) => AsBearing p f Double where
+  _Bearing =
+    prism'
+      (\(Bearing i) -> i)
+      (\i -> case i >= 0 && i < 360 of
+               True -> Just (Bearing i)
+               False -> Nothing)
+
diff --git a/src/Data/Geo/Geodetic/Curve.hs b/src/Data/Geo/Geodetic/Curve.hs
--- a/src/Data/Geo/Geodetic/Curve.hs
+++ b/src/Data/Geo/Geodetic/Curve.hs
@@ -1,16 +1,22 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GADTs #-}
+
 -- | A geodetic curve is made of a distance in metres, an azimuth and a reverse azimuth.
 module Data.Geo.Geodetic.Curve(
   Curve
+, AsCurve(..)
 , curve
 , curveDistance
 , curveAzimuth
 , curveReverseAzimuth
 ) where
 
-import Prelude(Eq, Show(..), Ord(..), Double, (.), abs, showString, showParen)
+import Prelude(Eq, Show(..), Ord(..), Double, showString, showParen, id)
+import Data.Functor(Functor)
 import Data.List(unwords)
 import Text.Printf(printf)
-import Control.Lens(Lens', lens)
+import Control.Lens(Optic', Lens', Profunctor, lens, iso)
 import Data.Geo.Geodetic.Azimuth
 
 data Curve =
@@ -49,3 +55,17 @@
   Lens' Curve Azimuth
 curveReverseAzimuth =
   lens (\(Curve _ _ r) -> r) (\(Curve d a _) r -> Curve d a r)
+
+class AsCurve p f s where
+  _Curve ::
+    Optic' p f s Curve
+
+instance AsCurve p f Curve where
+  _Curve =
+    id
+
+instance (Profunctor p, Functor f) => AsCurve p f (Double, Azimuth, Azimuth) where
+  _Curve =
+    iso
+      (\(d, a, r) -> Curve d a r)
+      (\(Curve d a r) -> (d, a, r))
diff --git a/src/Data/Geo/Geodetic/Ellipsoid.hs b/src/Data/Geo/Geodetic/Ellipsoid.hs
--- a/src/Data/Geo/Geodetic/Ellipsoid.hs
+++ b/src/Data/Geo/Geodetic/Ellipsoid.hs
@@ -1,12 +1,16 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GADTs #-}
+
 module Data.Geo.Geodetic.Ellipsoid(
   -- * Data type
   Ellipsoid
-, HasEllipsoid(..)
+, AsEllipsoid(..)
   -- * Ellipsoid properties
-, HasSemiMajor(..)
-, HasSemiMinor(..)
-, HasFlattening(..)
-, HasInverseFlattening(..)
+, AsSemiMajor(..)
+, AsSemiMinor(..)
+, AsFlattening(..)
+, AsInverseFlattening(..)
   -- * Ellipsoid construction
 , semiMajorFlattening
 , semiMinorFlattening
@@ -30,8 +34,8 @@
 , clarke1880
 ) where
 
-import Prelude(Eq, Ord, Show, Num(..), Fractional(..), Double, id)
-import Control.Lens(Lens', lens)
+import Prelude(Eq, Ord, Show, Num(..), Fractional(..), Double, id, Functor)
+import Control.Lens(Optic', lens)
 
 data Ellipsoid =
   Ellipsoid
@@ -67,61 +71,61 @@
 semiMinorInverseFlattening s f =
   semiMinorFlattening s (1/f)
 
-class HasEllipsoid t where
-  ellipsoid ::
-    Lens' t Ellipsoid
+class AsEllipsoid p f s where
+  _Ellipsoid ::
+    Optic' p f s Ellipsoid    
 
-instance HasEllipsoid Ellipsoid where
-  ellipsoid =
+instance AsEllipsoid p f Ellipsoid where
+  _Ellipsoid =
     id
 
-class HasSemiMajor t where
-  semiMajor ::
-    Lens' t Double
+class AsSemiMajor p f s where
+  _SemiMajor ::
+    Optic' p f s Double
 
-instance HasSemiMajor Double where
-  semiMajor =
+instance AsSemiMajor p f Double where
+  _SemiMajor =
     id
 
-instance HasSemiMajor Ellipsoid where
-  semiMajor =
-    lens (\(Ellipsoid s _) -> s) (\(Ellipsoid _ f) s -> Ellipsoid s f)
+instance (p ~ (->), Functor f) => AsSemiMajor p f Ellipsoid where
+  _SemiMajor =
+     lens (\(Ellipsoid s _) -> s) (\(Ellipsoid _ f) s -> Ellipsoid s f)   
 
-class HasSemiMinor t where
-  semiMinor ::
-    Lens' t Double
+class AsSemiMinor p f s where
+  _SemiMinor ::
+    Optic' p f s Double
 
-instance HasSemiMinor Double where
-  semiMinor =
-    id
+instance AsSemiMinor p f Double where
+  _SemiMinor =
+    id    
 
-instance HasSemiMinor Ellipsoid where
-  semiMinor =
-    lens (\(Ellipsoid s f) -> (1.0 - f) * s) (\(Ellipsoid _ f) s -> Ellipsoid (s * f - 1.0) f)
+instance (p ~ (->), Functor f) => AsSemiMinor p f Ellipsoid where
+  _SemiMinor =
+     lens (\(Ellipsoid s f) -> (1.0 - f) * s) (\(Ellipsoid _ f) s -> Ellipsoid (s * f - 1.0) f)
 
-class HasFlattening t where
-  flattening ::
-    Lens' t Double
+class AsFlattening p f s where
+  _Flattening ::
+    Optic' p f s Double
 
-instance HasFlattening Double where
-  flattening =
+instance AsFlattening p f Double where
+  _Flattening =
     id
 
-instance HasFlattening Ellipsoid where
-  flattening =
+instance (p ~ (->), Functor f) => AsFlattening p f Ellipsoid where
+  _Flattening =
     lens (\(Ellipsoid _ f) -> f) (\(Ellipsoid s _) f -> Ellipsoid s f)
 
-class HasInverseFlattening t where
-  inverseFlattening ::
-    Lens' t Double
+class AsInverseFlattening p f s where
+  _InverseFlattening ::
+    Optic' p f s Double
 
-instance HasInverseFlattening Double where
-  inverseFlattening =
+instance AsInverseFlattening p f Double where
+  _InverseFlattening =
     id
 
-instance HasInverseFlattening Ellipsoid where
-  inverseFlattening =
-    lens (\(Ellipsoid _ f) -> 1/f) (\(Ellipsoid s _) f -> Ellipsoid s (1/f))
+instance (p ~ (->), Functor f) => AsInverseFlattening p f Ellipsoid where
+  _InverseFlattening =
+    lens (\(Ellipsoid _ f) -> f) (\(Ellipsoid s _) f -> Ellipsoid s f)
 
 wgs84 ::
   Ellipsoid
diff --git a/src/Data/Geo/Geodetic/GreatCircle.hs b/src/Data/Geo/Geodetic/GreatCircle.hs
--- a/src/Data/Geo/Geodetic/GreatCircle.hs
+++ b/src/Data/Geo/Geodetic/GreatCircle.hs
@@ -8,6 +8,7 @@
 ) where
 
 import Prelude(Double, Num(..), Fractional(..), pi, sin, cos, acos)
+import Control.Applicative(Const)
 import Control.Lens((#), (^.))
 import System.Args.Optional(Optional1(..))
 import Data.Geo.Coordinate
@@ -20,54 +21,54 @@
 
 -- | Great circle spherical law algorithm.
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLaw earthMean fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (sphericalLaw earthMean fr to)) :: Maybe String
 -- Just "15000950.5589"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw earthMean fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (sphericalLaw earthMean fr to)) :: Maybe String
 -- Just "17128743.0669"
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLaw (6350000 ^. nSphere) fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (sphericalLaw ((6350000 :: Double) ^. _Sphere) fr to)) :: Maybe String
 -- Just "14959840.4461"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw (6350000 ^. nSphere) fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (sphericalLaw ((6350000 :: Double) ^. _Sphere) fr to)) :: Maybe String
 -- Just "17081801.7377"
 sphericalLaw ::
-  (HasCoordinate c1, HasCoordinate c2) =>
+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>
   Sphere -- ^ reference sphere
-  -> c1 -- ^ start coordinate
-  -> c2 -- ^ end coordinate
+  -> start -- ^ start coordinate
+  -> end -- ^ end coordinate
   -> Double
 sphericalLaw s start' end' =
-  let start = start' ^. coordinate
-      end = end' ^. coordinate
+  let start = start' ^. _Coordinate
+      end = end' ^. _Coordinate
       toRadians n = n * pi / 180
-      lat1 = toRadians (fracLatitude # (start ^. latitude))
-      lat2 = toRadians (fracLatitude # (end ^. latitude))
-      lon1 = toRadians (fracLongitude # (start ^. longitude))
-      lon2 = toRadians (fracLongitude # (end ^. longitude))
-  in acos (sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)) * nSphere # s
+      lat1 = toRadians (_Latitude # (start ^. _Latitude))
+      lat2 = toRadians (_Latitude # (end ^. _Latitude))
+      lon1 = toRadians (_Longitude # (start ^. _Longitude))
+      lon2 = toRadians (_Longitude # (end ^. _Longitude))
+  in acos (sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)) * _Sphere # s
 
 -- | Great circle spherical law algorithm with a default sphere of the earth mean.
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLawD fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (sphericalLawD fr to)) :: Maybe String
 -- Just "15000950.5589"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLawD fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (sphericalLawD fr to)) :: Maybe String
 -- Just "17128743.0669"
 sphericalLawD ::
-  (HasCoordinate c1, HasCoordinate c2) =>
-  c1 -- ^ start coordinate
-  -> c2 -- ^ end coordinate
+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>
+  start -- ^ start coordinate
+  -> end -- ^ end coordinate
   -> Double
 sphericalLawD =
   sphericalLaw earthMean
 
 -- | Great circle spherical law algorithm with an optionally applied default sphere of the earth mean.
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLaw' fr to :: Double)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (sphericalLaw' fr to :: Double)) :: Maybe String
 -- Just "15000950.5589"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw' fr to :: Double)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (sphericalLaw' fr to :: Double)) :: Maybe String
 -- Just "17128743.0669"
 sphericalLaw' ::
   (Optional1
diff --git a/src/Data/Geo/Geodetic/Haversine.hs b/src/Data/Geo/Geodetic/Haversine.hs
--- a/src/Data/Geo/Geodetic/Haversine.hs
+++ b/src/Data/Geo/Geodetic/Haversine.hs
@@ -1,4 +1,6 @@
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 -- | Haversine geodetic distance algorithm.
 module Data.Geo.Geodetic.Haversine(
@@ -7,71 +9,72 @@
 , haversine'
 ) where
 
-import Prelude(Double, Num(..), Fractional(..), (.), pi, atan, sin, atan2, cos, sqrt)
+import Control.Applicative(Const)
+import Prelude(Double, Num(..), Fractional(..), (.), pi, sin, atan2, cos, sqrt)
 import Control.Lens((#), (^.))
 import System.Args.Optional(Optional1(..))
 import Data.Geo.Coordinate
 import Data.Geo.Geodetic.Sphere
 
 -- $setup
--- >>> import Prelude(Functor(..), Monad(..), String)
+-- >>> import Prelude(Functor(..), Monad(..), String, Double)
 -- >>> import Data.Maybe(Maybe)
 -- >>> import Text.Printf(printf)
 
 -- | Haversine algorithm.
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversine earthMean fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (haversine earthMean fr to)) :: Maybe String
 -- Just "15000950.5589"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine earthMean fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (haversine earthMean fr to)) :: Maybe String
 -- Just "17128743.0669"
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversine (6350000 ^. nSphere) fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (haversine ((6350000 :: Double) ^. _Sphere) fr to)) :: Maybe String
 -- Just "14959840.4461"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine (6350000 ^. nSphere) fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (haversine ((6350000 :: Double) ^. _Sphere) fr to)) :: Maybe String
 -- Just "17081801.7377"
 haversine ::
-  (HasCoordinate c1, HasCoordinate c2) =>
+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>
   Sphere -- ^ reference sphere
-  -> c1 -- ^ start coordinate
-  -> c2 -- ^ end coordinate
+  -> start -- ^ start coordinate
+  -> end -- ^ end coordinate
   -> Double
 haversine s start' end' =
-  let start = start' ^. coordinate
-      end = end' ^. coordinate
-      lat1 = fracLatitude # (start ^. latitude)
-      lat2 = fracLatitude # (end ^. latitude)
+  let start = start' ^. _Coordinate
+      end = end' ^. _Coordinate
+      lat1 = _Latitude # (start ^. _Latitude)
+      lat2 = _Latitude # (end ^. _Latitude)
       toRadians n = n * pi / 180
       dlat = (toRadians (lat1 - lat2)) / 2
-      dlon = (toRadians (fracLongitude # (start ^. longitude) - fracLongitude # (end ^. longitude))) / 2
+      dlon = (toRadians (_Longitude # (start ^. _Longitude) - _Longitude # (end ^. _Longitude))) / 2
       cosr = cos . toRadians
       square x = x * x
       a = square (sin dlat) + cosr lat1 * cosr lat2 * square (sin (dlon))
       c = 2 * atan2 (sqrt a) (sqrt (1 - a))
-  in (nSphere # s) * c
+  in (_Sphere # s) * c
 
 -- | Haversine algorithm with a default sphere of the earth mean.
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversineD fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (haversineD fr to)) :: Maybe String
 -- Just "15000950.5589"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversineD fr to)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (haversineD fr to)) :: Maybe String
 -- Just "17128743.0669"
 haversineD ::
-  (HasCoordinate c1, HasCoordinate c2) =>
-  c1 -- ^ start coordinate
-  -> c2 -- ^ end coordinate
+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>
+  start -- ^ start coordinate
+  -> end -- ^ end coordinate
   -> Double
 haversineD =
   haversine earthMean
 
 -- | Haversine algorithm with an optionally applied default sphere of the earth mean.
 --
--- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversine' fr to :: Double)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (haversine' fr to :: Double)) :: Maybe String
 -- Just "15000950.5589"
 --
--- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine' fr to :: Double)) :: Maybe String
+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (haversine' fr to :: Double)) :: Maybe String
 -- Just "17128743.0669"
 haversine' ::
   (Optional1
diff --git a/src/Data/Geo/Geodetic/Sphere.hs b/src/Data/Geo/Geodetic/Sphere.hs
--- a/src/Data/Geo/Geodetic/Sphere.hs
+++ b/src/Data/Geo/Geodetic/Sphere.hs
@@ -1,13 +1,17 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+
+
 -- | A sphere with a radius in metres.
 module Data.Geo.Geodetic.Sphere(
   Sphere
-, HasSphere(..)
-, nSphere
+, AsSphere(..)
 , earthMean
 ) where
 
-import Prelude(Double, Bool(..), Eq, Show(..), Ord(..), id, (&&), (++), (.), abs, showParen, showString)
-import Control.Lens(Iso', Lens', iso)
+import Prelude(Double, Eq, Show(..), Ord(..), id, (++), showParen, showString)
+import Control.Lens(Optic', Profunctor, iso)
+import Data.Functor(Functor)
 import Text.Printf(printf)
 
 -- $setup
@@ -25,30 +29,31 @@
   showsPrec n (Sphere d) =
     showParen (n > 10) (showString ("Sphere " ++ printf "%0.4f" d))
 
--- | An isomorphism on sphere to a double.
---
--- >>> 7 ^. nSphere
--- Sphere 7.0000
---
--- >>> 0 ^. nSphere
--- Sphere 0.0000
---
--- >>> (-7) ^. nSphere
--- Sphere -7.0000
-nSphere ::
-  Iso' Double Sphere
-nSphere =
-  iso Sphere (\(Sphere d) -> d)
-
 earthMean ::
   Sphere
 earthMean =
   Sphere 6367450
 
-class HasSphere t where
-  sphere ::
-    Lens' t Sphere
+class AsSphere p f s where
+  _Sphere ::
+    Optic' p f s Sphere
 
-instance HasSphere Sphere where
-  sphere =
+instance AsSphere p f Sphere where
+  _Sphere =
     id
+
+-- | An isomorphism on sphere to a double.
+--
+-- >>> (7 :: Double) ^. _Sphere
+-- Sphere 7.0000
+--
+-- >>> (0 :: Double) ^. _Sphere
+-- Sphere 0.0000
+--
+-- >>> (-7 :: Double) ^. _Sphere
+-- Sphere -7.0000
+instance (Functor f, Profunctor p) => AsSphere p f Double where
+  _Sphere =
+    iso 
+      Sphere
+      (\(Sphere d) -> d)
diff --git a/src/Data/Geo/Geodetic/Vincenty.hs b/src/Data/Geo/Geodetic/Vincenty.hs
--- a/src/Data/Geo/Geodetic/Vincenty.hs
+++ b/src/Data/Geo/Geodetic/Vincenty.hs
@@ -1,4 +1,7 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE GADTs #-}
 
 -- | An implementation of Thaddeus Vincenty's direct and inverse geodetic algorithms. <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>
 module Data.Geo.Geodetic.Vincenty(
@@ -13,8 +16,10 @@
 , inverse'
 ) where
 
-import Prelude(Eq(..), Show(..), Ord(..), Num(..), Floating(..), Fractional(..), Double, Int, Bool, Ordering(..), subtract, cos, sin, asin, tan, sqrt, atan, atan2, pi, (.), (++), (&&), ($!), error)
-import Control.Lens(lens, (^.), (#), (^?))
+import Prelude(Eq(..), Show(..), Ord(..), Num(..), Floating(..), Fractional(..), Double, Int, Bool, Ordering(..), subtract, cos, sin, asin, tan, sqrt, atan, atan2, pi, (.), (++), (&&), ($!), error, id)
+import Control.Applicative(Const)
+import Control.Lens(Profunctor, Prism', Optic', Iso', (^.), (#), (^?), iso, _1, _2, from)
+import Data.Functor(Functor)
 import Data.Maybe(fromMaybe)
 import System.Args.Optional(Optional2(..))
 import Data.Geo.Coordinate
@@ -41,29 +46,43 @@
     Bearing
   deriving (Eq, Ord, Show)
 
-instance HasCoordinate VincentyDirectResult where
-  coordinate =
-    lens (\(VincentyDirectResult c _) -> c) (\(VincentyDirectResult _ b) c -> VincentyDirectResult c b)
+class AsVincentyDirectResult p f s where
+  _VincentyDirectResult ::
+    Optic' p f s VincentyDirectResult
 
-instance HasBearing VincentyDirectResult where
-  bearing =
-    lens (\(VincentyDirectResult _ b) -> b) (\(VincentyDirectResult c _) b -> VincentyDirectResult c b)
+instance AsVincentyDirectResult p f VincentyDirectResult where
+  _VincentyDirectResult =
+    id
 
+instance (Profunctor p, Functor f) => AsVincentyDirectResult p f (Coordinate, Bearing) where
+  _VincentyDirectResult =
+    iso
+      (\(c, b) -> VincentyDirectResult c b)
+      (\(VincentyDirectResult c b) -> (c, b))
+
+instance (p ~ (->), Functor f) => AsCoordinate p f VincentyDirectResult where
+  _Coordinate =
+    from (_VincentyDirectResult :: Iso' (Coordinate, Bearing) VincentyDirectResult) . _1
+
+instance (p ~ (->), Functor f) => AsBearing p f VincentyDirectResult where
+  _Bearing =
+    from (_VincentyDirectResult :: Iso' (Coordinate, Bearing) VincentyDirectResult) . _2
+
 -- | Vincenty direct algorithm.
 --
--- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)
+-- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) (27.812 <°> 154.295)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0981)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1466))) (Bearing 165.3451))
 --
--- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84)
+-- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) ((-66.093) <°> 12.84)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183))
 --
--- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)
+-- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) (27.812 <°> 154.295)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0986)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1464))) (Bearing 165.3451))
 --
--- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84)
+-- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) ((-66.093) <°> 12.84)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0662)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4139))) (Bearing 165.3183))
 direct ::
-  (HasEllipsoid e, HasCoordinate c, HasBearing b) =>
+  (AsCoordinate (->) (Const Coordinate) c, AsBearing (->) (Const Bearing) b, AsEllipsoid (->) (Const Ellipsoid) e) =>
   e -- ^ reference ellipsoid
   -> Convergence -- ^ convergence point to stop calculating
   -> c -- ^ begin coordinate
@@ -71,22 +90,24 @@
   -> Double -- ^ distance
   -> VincentyDirectResult
 direct e' conv start' bear' dist =
-  let e = e' ^. ellipsoid
-      start = start' ^. coordinate
-      bear = bear' ^. bearing
-      sMnr = e ^. semiMinor
-      flat = e ^. flattening
+  let radianLatitude :: Prism' Double Latitude
+      radianLatitude = iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Latitude
+      e = e' ^. _Ellipsoid
+      start = start' ^. _Coordinate
+      bear = bear' ^. _Bearing
+      sMnr = e ^. _SemiMinor
+      flat = e ^. _Flattening
       alpha = radianBearing # bear
       cosAlpha = cos alpha
       sinAlpha = sin alpha
-      tanu1 = (1.0 - flat) * tan (radianLatitude # (start ^. latitude))
+      tanu1 = (1.0 - flat) * tan (radianLatitude # (start ^. _Latitude))
       cosu1 = 1.0 / sqrt (1.0 + square tanu1)
       sinu1 = tanu1 * cosu1
       sigma1 = atan2 tanu1 cosAlpha
       csa = cosu1 * sinAlpha
       sin2Alpha = square csa
       cos2Alpha = 1 - sin2Alpha
-      ab d f g h i = let s = cos2Alpha * (square (e ^. semiMajor / sMnr) - 1)
+      ab d f g h i = let s = cos2Alpha * (square (e ^. _SemiMajor / sMnr) - 1)
                      in (s / d) * (f + s * (g + s * (h - i * s)))
       a = 1 + ab 16384 4096 (-768) 320 175
       b = ab 1024 256 (-128) 74 47
@@ -111,8 +132,8 @@
       sss = sinu1 * sinSigma
       latitude' = let r = atan2 (sinu1 * cosSigma + cosu1 * sinSigma * cosAlpha) ((1.0 - flat) * sqrt (sin2Alpha + (sss - ccca) ** 2.0))
                   in fromMaybe (error ("Invariant not met. Latitude in radians not within range " ++ show r)) (r ^? radianLatitude)
-      longitude' = let r = fracLongitude # (start ^. longitude) + ((atan2 (sinSigma * sinAlpha) (cc - sss * cosAlpha) - (1 - c) * flat * csa * (sigma'' + c * sinSigma * (cosSigmaM2 + c * cosSigma * (-1 + 2 * cos2SigmaM2)))) * 180 / pi)
-                    in fromMaybe (error ("Invariant not met. Longitude in radians not within range " ++ show r)) (r ^? fracLongitude)
+      longitude' = let r = _Longitude # (start ^. _Longitude) + ((atan2 (sinSigma * sinAlpha) (cc - sss * cosAlpha) - (1 - c) * flat * csa * (sigma'' + c * sinSigma * (cosSigmaM2 + c * cosSigma * (-1 + 2 * cos2SigmaM2)))) * 180 / pi)
+                   in fromMaybe (error ("Invariant not met. Longitude in radians not within range " ++ show r)) (r ^? _Longitude)
   in VincentyDirectResult
        (latitude' .#. longitude')
        (
@@ -122,13 +143,13 @@
 
 -- | Vincenty direct algorithm with a default ellipsoid of WGS84 and standard convergence.
 --
--- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)
+-- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) (27.812 <°> 154.295)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0981)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1466))) (Bearing 165.3451))
 --
--- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84)
+-- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) ((-66.093) <°> 12.84)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183))
 directD ::
-  (HasCoordinate c, HasBearing b) =>
+  (AsCoordinate (->) (Const Coordinate) c, AsBearing (->) (Const Bearing) b) =>
   c -- ^ begin coordinate
   -> b -- ^ bearing
   -> Double -- ^ distance
@@ -138,10 +159,10 @@
 
 -- | Vincenty direct algorithm with an optionally applied default ellipsoid of WGS84 and standard convergence.
 --
--- >>> fmap (\c' -> direct' c' (modBearing 165.34) (4235 :: Double) :: VincentyDirectResult) (27.812 ..#.. 154.295)
+-- >>> fmap (\c' -> direct' c' (modBearing 165.34) (4235 :: Double) :: VincentyDirectResult) (27.812 <°> 154.295)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0981)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1466))) (Bearing 165.3451))
 --
--- >>> fmap (\c' -> direct' c' (modBearing 165.34) (4235 :: Double) :: VincentyDirectResult) ((-66.093) ..#.. 12.84)
+-- >>> fmap (\c' -> direct' c' (modBearing 165.34) (4235 :: Double) :: VincentyDirectResult) ((-66.093) <°> 12.84)
 -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183))
 direct' ::
   (Optional2
@@ -159,38 +180,42 @@
 
 -- | Vincenty inverse algorithm.
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse wgs84 convergence fr to)
+-- >>> do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (inverse wgs84 convergence fr to)
 -- Just (GeodeticCurve 14998576.9860 Azimuth 180.0000 Azimuth 0.0000)
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse wgs84 convergence fr to)
+-- >>> do fr <- 27.812 <°> 154.295; to <- 87.7769 <°> 19.944; return (inverse wgs84 convergence fr to)
 -- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000)
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse ans convergence fr to)
+-- >>> do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (inverse ans convergence fr to)
 -- Just (GeodeticCurve 14998630.4056 Azimuth 180.0000 Azimuth 0.0000)
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse ans convergence fr to)
+-- >>> do fr <- 27.812 <°> 154.295; to <- 87.7769 <°> 19.944; return (inverse ans convergence fr to)
 -- Just (GeodeticCurve 7099229.9126 Azimuth 0.0000 Azimuth 180.0000)
 inverse ::
-  (HasEllipsoid e, HasCoordinate c1, HasCoordinate c2) =>
+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end, AsEllipsoid (->) (Const Ellipsoid) e) =>
   e -- ^ reference ellipsoid
   -> Convergence -- ^ convergence point to stop calculating
-  -> c1 -- ^ start coordinate
-  -> c2 -- ^ end coordinate
+  -> start -- ^ start coordinate
+  -> end -- ^ end coordinate
   -> Curve
 inverse e' conv start' end' =
-  let e = e' ^. ellipsoid
-      start = start' ^. coordinate
-      end = end' ^. coordinate
-      b = e ^. semiMinor
-      f = e ^. flattening
+  let radianLatitude :: Prism' Double Latitude
+      radianLatitude = iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Latitude
+      radianLongitude :: Prism' Double Longitude
+      radianLongitude = iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Longitude
+      e = e' ^. _Ellipsoid
+      start = start' ^. _Coordinate
+      end = end' ^. _Coordinate
+      b = e ^. _SemiMinor
+      f = e ^. _Flattening
       (phi1, phi2) =
-        let rl k = radianLatitude # (k ^. latitude)
+        let rl k = radianLatitude # (k ^. _Latitude)
         in (rl start, rl end)
       a2b2b2 =
         let ss z = square (z e)
-        in ss (^. semiMajor) / ss (^. semiMinor) - 1
+        in ss (^. _SemiMajor) / ss (^. _SemiMinor) - 1
       omega =
-        let rl k = radianLongitude # (k ^. longitude)
+        let rl k = radianLongitude # (k ^. _Longitude)
         in rl end - rl start
       (u1, u2) =
         let at = atan . ((1 - f) *) . tan
@@ -246,25 +271,25 @@
 
 -- | Vincenty inverse algorithm with a default ellipsoid of WGS84 and standard convergence.
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverseD fr to)
+-- >>> do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (inverseD fr to)
 -- Just (GeodeticCurve 14998576.9860 Azimuth 180.0000 Azimuth 0.0000)
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverseD fr to)
+-- >>> do fr <- 27.812 <°> 154.295; to <- 87.7769 <°> 19.944; return (inverseD fr to)
 -- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000)
 inverseD ::
-  (HasCoordinate c1, HasCoordinate c2) =>
-  c1 -- ^ start coordinate
-  -> c2 -- ^ end coordinate
+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>
+  start -- ^ start coordinate
+  -> end -- ^ end coordinate
   -> Curve
 inverseD =
   inverse wgs84 convergence
 
 -- | Vincenty inverse algorithm with an optionally applied default ellipsoid of WGS84 and standard convergence.
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse' fr to :: Curve)
+-- >>> do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (inverse' fr to :: Curve)
 -- Just (GeodeticCurve 14998576.9860 Azimuth 180.0000 Azimuth 0.0000)
 --
--- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse' fr to :: Curve)
+-- >>> do fr <- 27.812 <°> 154.295; to <- 87.7769 <°> 19.944; return (inverse' fr to :: Curve)
 -- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000)
 inverse' ::
   (Optional2