2023-10-09 16:12:54 (UTC-03:00)

Marcel Rodrigues <marcelgmr@gmail.com>

map: add Mercator projection

diff --git a/lib/anim/map.lua b/lib/anim/map.lua
index b1fa50e..048b442 100644
--- a/lib/anim/map.lua
+++ b/lib/anim/map.lua
@@ -1,11 +1,12 @@
 local util = require "anim.util"
 local bio = require "anim.bio"
 
-local abs = math.abs
+local abs, min, max = math.abs, math.min, math.max
 local rad, deg = math.rad, math.deg
 local cos, sin, tan = math.cos, math.sin, math.tan
+local cosh, sinh, tanh = math.cosh, math.sinh, math.tanh
 local acos, asin, atan, atan2 = math.acos, math.asin, math.atan, math.atan2
-local sqrt = math.sqrt
+local sqrt, log = math.sqrt, math.log
 local pi, huge = math.pi, math.huge
 
 -- == Utilities ==
@@ -94,8 +95,43 @@ local function centroid(region)
     return lon1, lat1
 end
 
+-- correct angle such that -pi <= a <= pi
+local function corangle(a)
+    while a < -pi do
+        a = a + 2 * pi
+    end
+    while a > pi do
+        a = a - 2 * pi
+    end
+    return a
+end
+
 -- == Projections ==
 
+-- Mercator Projection for the Spherical Earth.
+
+local Mercator = {}
+Mercator.__index = Mercator
+
+function Mercator:map(lon, lat)
+    lat = max(lat, -85.5) -- avoid distortions near singularity at latitude -90
+    lon, lat = rad(lon), rad(lat)
+    local x = self.r * corangle(lon - self.lon)
+    local y = self.r * log(tan(pi/4 + lat/2))
+    return x, y
+end
+
+function Mercator:inv(x, y)
+    local lon = corangle(x / self.r + self.lon)
+    local lat = atan(sinh(y / self.r))
+    lon, lat = deg(lon), deg(lat)
+    return lon, lat
+end
+
+function Mercator:model()
+    return {type="sphere", r=self.r}
+end
+
 -- Lambert Azimuthal Equal-Area Projection for the Spherical Earth.
 
 local AzimuthalEqualArea = {}
@@ -137,7 +173,8 @@ function AzimuthalEqualArea:model()
 end
 
 local projs = {
-    AzimuthalEqualArea=AzimuthalEqualArea
+    Mercator=Mercator,
+    AzimuthalEqualArea=AzimuthalEqualArea,
 }
 
 -- Generic projection interface.