spherical (empty) → 0.1.0.0
raw patch · 6 files changed
+249/−0 lines, 6 filesdep +basedep +composition-prelude
Dependencies added: base, composition-prelude
Files
- CHANGELOG.md +5/−0
- LICENSE +11/−0
- README.md +4/−0
- cabal.project.local +2/−0
- spherical.cabal +47/−0
- src/Math/Geometry/Spherical.hs +180/−0
+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# spherical++## 0.1.0.0++Initial release
+ LICENSE view
@@ -0,0 +1,11 @@+Copyright Vanessa McHale (c) 2018++Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.++3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,4 @@+# spherical++A library for working with geometry on the surface of the sphere, mainly for use+in GIS.
+ cabal.project.local view
@@ -0,0 +1,2 @@+constraints: spherical +development+documentation: true
+ spherical.cabal view
@@ -0,0 +1,47 @@+cabal-version: 1.18+name: spherical+version: 0.1.0.0+license: BSD3+license-file: LICENSE+copyright: Copyright: (c) 2018 Vanessa McHale+maintainer: vanessa.mchale@iohk.io+author: Vanessa McHale+synopsis: Geometry on a sphere+description:+ Functions in pure Haskell for geometric computations on the surface of a sphere.+category: Math+build-type: Simple+extra-source-files:+ cabal.project.local+extra-doc-files: README.md+ CHANGELOG.md++source-repository head+ type: git+ location: https://github.com/vmchale/spherical++flag development+ description:+ Enable `-Werror`+ default: False+ manual: True++library+ exposed-modules:+ Math.Geometry.Spherical+ hs-source-dirs: src+ default-language: Haskell98+ ghc-options: -Wall+ build-depends:+ base >=4.3 && <5,+ composition-prelude -any+ + if flag(development)+ ghc-options: -Werror+ + if impl(ghc >=8.0)+ ghc-options: -Wincomplete-uni-patterns -Wincomplete-record-updates+ -Wredundant-constraints -Widentities+ + if impl(ghc >=8.4)+ ghc-options: -Wmissing-export-lists
+ src/Math/Geometry/Spherical.hs view
@@ -0,0 +1,180 @@+{-# LANGUAGE RankNTypes #-}++module Math.Geometry.Spherical+ ( -- * Computations+ areaTriangle+ , compactness+ , relativeCompactness+ , areaConvex+ , perimeterPolygon+ , areaPolygon+ -- * Projections+ , albers+ , littow+ , craig+ , winkel3+ , mercator+ -- * Reference points+ , washingtonDC+ , mecca+ -- * Helper functions+ , radians+ , toRadians+ , project+ ) where++import Control.Composition++-- | Convert from degrees to radians.+radians :: Floating a => a -> a+radians = (*(pi/180))++-- | Convert both coördinates to radians.+toRadians :: Floating a => (a, a) -> (a, a)+toRadians = both radians++sinc :: Floating a => a -> a+sinc x = sin x / x++-- | For use as a reference point+washingtonDC :: (Floating a) => (a, a)+washingtonDC = toRadians (38.9072, -77.0369)++-- | For use as a reference point+mecca :: (Floating a) => (a, a)+mecca = toRadians (21.3891, 39.8579)++-- | Littow retroazimuthal + conformal projection+littow :: Floating a => (a, a) -> (a, a)+littow (long, lat) = (sin (long - referenceLong) /cos lat, cos (long - referenceLong) * tan lat)+ where referenceLong = radians 0++-- | Craig retroazimuthal projection+craig :: (Floating a, Eq a)+ => (a, a) -- ^ Reference point given in radians+ -> ((a, a) -> (a, a))+craig referencePoint (long, lat) = (long - referenceLong, y)+ where (referenceLong, referenceLat) = referencePoint+ expr = sin lat * cos (long - referenceLong) - tan referenceLat * cos lat+ y | long - referenceLong == 0 = expr+ | otherwise = (long - referenceLong) / sin (long - referenceLong) * expr++-- | Winkel Tripel projection+winkel3 :: Floating a => (a, a) -> (a, a)+winkel3 (long, lat) = ((lambda * cos phi1 + (2 * cos lat * sin (lambda/2)/sinc alpha))/2, (lat + sin lat/sinc alpha)/2)+ where lambda = long - lambda0+ phi1 = acos $ 2 / pi+ alpha = acos $ cos lat * cos (lambda/2)+ lambda0 = radians (-77.0369)++-- | Mercator projection.+mercator :: Floating a => (a, a) -> (a, a)+mercator (long, lat) = (long - meridian, asinh (tan lat))+ where meridian = radians (-98.5795)++-- | Bonne projection.+bonne :: Floating a+ => a -- ^ Standard Parallel. If you are unsure of what to put, try @radians 45@.+ -> a -- ^ Central meridian. If you are unsure of what to put, try @snd washingtonDC@.+ -> (a, a)+ -> (a, a)+bonne phi1 meridian (long, lat) = (rho * sin e, cot phi1 - rho * cos e)+ where rho = cot phi1 + phi1 - lat+ e = (long - meridian) * cos lat / rho+ cot = (1/) . tan++-- | Albers projection for a given reference point.+--+-- > ablers washingtonDC+albers :: Floating a => (a, a) -- ^ A reference point on the sphere+ -> ((a, a) -> (a, a))+albers referencePoint (long, lat) = (rho * sin theta, rho' - rho * cos theta)+ where n = (sin phi1 + sin phi2) / 2+ theta = n * (long - referenceLong)+ c = cos phi1^(2 :: Int) + 2 * n * sin phi1+ rho = sqrt (c - 2 * n * sin lat) / n+ rho' = sqrt (c - 2 * n * sin referenceLat) / n+ -- standard parallels @ 20, 50 degrees+ phi1 = radians 20+ phi2 = radians 50+ (referenceLong, referenceLat) = referencePoint++-- | Project a given polygon+project :: (Floating a, Functor f) => ((a, a) -> (a, a)) -- ^ A projection+ -> f (a, a) -- ^ A polygon defined by points on the sphere, in degrees+ -> f (a, a)+project f = fmap (f . toRadians)++-- | Compute the area of a triangle using L\'Huillier\'s formula+areaTriangle :: Floating a => a -- ^ Radius of the sphere+ -> (a, a) -- ^ A point given in radians+ -> (a, a)+ -> (a, a)+ -> a+areaTriangle r x1 x2 x3 = r^(2 :: Int) * e+ where e = 4 * atan(sqrt(tan(s/2) * tan((s - a)/2) * tan((s - b)/2) * tan((s - c)/2)))+ s = (a + b + c) / 2+ a = distanceRad x1 x2+ b = distanceRad x1 x3+ c = distanceRad x2 x3+ distanceRad = on centralAngle toRadians++-- TODO mandelbrot/fractal dimension?+-- consider "area of largest circumscribable circle" as well.++-- | Relative compactness. Dimensionless.+relativeCompactness :: (Floating a) => [(a, a)] -> a+relativeCompactness = (*scale) . compactness+ where scale = 1/4*pi++-- | Take the area of the polygon and divide by the perimeter squared. Dimensionless.+compactness :: (Floating a)+ => [(a, a)] -- ^ Polygons on surface of the sphere, given in degrees+ -> a+compactness p = areaPolygon 1 p/(perimeterPolygon 1 p^(2 :: Int))++-- | Compute the area of a convex polygon on the surface of a sphere.+areaConvex :: (Floating a)+ => a -- ^ Radius of a sphere+ -> [(a, a)] -- ^ Polygon on the surface of the sphere+ -> a+areaConvex r (base1:base2:pts) = fst $ foldr stepArea (0,base2) pts+ where stepArea point (sum', base) = (sum' + areaTriangle r base1 base point, point)+areaConvex _ _ = error "attempted to take area of polygon with < 3 points"++-- | Uses areal projection; then finds area of the polygon by the shoelace+-- method.+areaPolygon :: Floating a+ => a -- ^ Radius of sphere+ -> [(a, a)] -- ^ Polygon on the sphere, with points given in degrees.+ -> a+areaPolygon r = (*factor) . areaPolyRectangular . fmap (bonne (radians 25) (snd washingtonDC) . toRadians)+ where factor = 1717856/4.219690791828533e-2 * ((r / 6371) ^ (2 :: Int))++perimeterPolygon :: Floating a+ => a -- ^ Radius of sphere+ -> [(a, a)] -- ^ Polygon on sphere given in degrees+ -> a+perimeterPolygon r [x1, x2] = distance r x1 x2+perimeterPolygon r (x1:x2:points) = perimeterPolygon r (x2:points) + distance r x1 x2+perimeterPolygon _ _ = error "Attempted to take perimeter of polygon with no points"++-- | Find the area of a polygon with rectangular coördinates given. This+-- function will throw an error when a polygon with no points is given.+areaPolyRectangular :: (Floating a) => [(a, a)] -> a+areaPolyRectangular (pt:pts) = abs . (*0.5) . fst $ foldr areaPolyCalc (0,pt) pts+ where areaPolyCalc (x2, y2) (sum',(x1,y1)) = (sum' + (x1 * y2 - x2 * y1),(x2,y2))+areaPolyRectangular _ = error "Attempted to take area of polygon with no points"++-- | Distance in kilometers between two points given in degrees.+distance :: Floating a+ => a -- ^ Radius of sphere+ -> (a, a) -- ^ Point on sphere given in degrees+ -> (a, a) -- ^ Point on sphere given in degrees+ -> a+distance r = (*r) .* on centralAngle toRadians++-- | Compute central angle from points given in radians+centralAngle :: Floating a => (a, a) -> (a, a) -> a+centralAngle (long1, lat1) (long2, lat2) =+ acos $ sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (long1 - long2)