2021-08-22 17:05:09 (UTC-03:00)

Marcel Rodrigues <marcelgmr@gmail.com>

add mapping utilities

diff --git a/map.lua b/map.lua
new file mode 100644
index 0000000..5e56549
--- /dev/null
+++ b/map.lua
@@ -0,0 +1,257 @@
+local abs = math.abs
+local rad, deg = math.rad, math.deg
+local cos, sin, tan = math.cos, math.sin, math.tan
+local acos, asin, atan, atan2 = math.acos, math.asin, math.atan, math.atan2
+local sqrt = math.sqrt
+local pi, huge = math.pi, math.huge
+
+-- == Utilities ==
+
+-- region is a list of polygons in geographic coordinates.
+
+local function distance(lon1, lat1, lon2, lat2, r)
+    r = r or 6378137
+    local dlat = rad(lat2 - lat1)
+    local dlon = rad(lon2 - lon1)
+    lat1, lat2 = rad(lat1), rad(lat2)
+    local a1, a2, a, c
+    a1 = sin(dlat/2) * sin(dlat/2)
+    a2 = sin(dlon/2) * sin(dlon/2) * cos(lat1) * cos(lat2)
+    a = a1 + a2
+    c = 2 * atan2(sqrt(a), sqrt(1-a))
+    return r * c
+end
+
+local function bbox(region)
+    local x0, y0, x1, y1 = huge, huge, -huge, -huge
+    for _, polygon in ipairs(region) do
+        for _, point in ipairs(polygon) do
+            local x, y = unpack(point)
+            x0 = x < x0 and x or x0
+            y0 = y < y0 and y or y0
+            x1 = x > x1 and x or x1
+            y1 = y > y1 and y or y1
+        end
+    end
+    return x0, y0, x1, y1
+end
+
+local function centroid(region)
+    local epsilon = 1e-10
+    local x0, y0, x1, y1 = bbox(region)
+    local lon0 = (x0 + x1) / 2
+    local lat0 = (y0 + y1) / 2
+    local lon1, lat1
+    while true do
+        local prj = Proj("AzimuthalEqualArea", {lon0, lat0})
+        local cw = {}
+        for i, polygon in ipairs(region) do
+            local xys = {}
+            for j, point in ipairs(polygon) do
+                xys[j] = {prj:map(unpack(point))}
+            end
+            if xys[#xys][0] ~= xys[1][0] or xys[#xys][1] ~= xys[1][1] then
+                xys[#xys+1] = xys[1]
+            end
+            -- http://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
+            local cx, cy, sa = 0, 0, 0
+            for j = 1, #xys-1 do
+                local x0, y0 = unpack(xys[j])
+                local x1, y1 = unpack(xys[j+1])
+                local f = x0 * y1 - x1 * y0
+                cx = cx + (x0 + x1) * f
+                cy = cy + (y0 + y1) * f
+                sa = sa + f
+            end
+            cx = cx / (3 * sa)
+            cy = cy / (3 * sa)
+            cw[#cw+1] = {cx, cy, sa}
+        end
+        local cx, cy, sw = 0, 0, 0
+        for i = 1, #cw do
+            local x, y, w = unpack(cw[i])
+            cx = cx + x * w
+            cy = cy + y * w
+            sw = sw + w
+        end
+        cx = cx / sw
+        cy = cy / sw
+        lon1, lat1 = prj:inv(cx, cy)
+        if abs(lon1-lon0) <= epsilon and abs(lat1-lat0) <= epsilon then
+            break
+        end
+        lon0, lat0 = lon1, lat1
+    end
+    return lon1, lat1
+end
+
+-- == Projections ==
+
+-- Lambert Azimuthal Equal-Area Projection for the Spherical Earth.
+
+local AzimuthalEqualArea = {}
+AzimuthalEqualArea.__index = AzimuthalEqualArea
+
+function AzimuthalEqualArea:map(lon, lat)
+    lon, lat = rad(lon), rad(lat)
+    lon = lon - self.lon0
+    local k, x, y
+    k = sqrt(2 / (1 + sin(self.lat0) * sin(lat) + cos(self.lat0) * cos(lat) * cos(lon)))
+    x = self.r * k * cos(lat) * sin(lon)
+    y = self.r * k * (cos(self.lat0) * sin(lat) - sin(self.lat0) * cos(lat) * cos(lon))
+    return x, y
+end
+
+function AzimuthalEqualArea:inv(x, y)
+    local p = sqrt(x*x + y*y)
+    local c = 2 * asin(p / (2 * self.r))
+    local lon, lat
+    -- FIXME: In the formulas below, should it be atan or atan2?
+    if self.lat0 == pi / 2 then
+        -- North Polar Aspect.
+        lon = self.lon0 + atan(x/(-y))
+    elseif self.lat0 == -pi / 2 then
+        -- South Polar Aspect.
+        lon = self.lon0 + atan(x/y)
+    else
+        -- Any other Oblique Aspect.
+        local den = p * cos(self.lat0) * cos(c) - y * sin(self.lat0) * sin(c)
+        lon = self.lon0 + atan(x * sin(c) / den)
+    end
+    lat = asin(cos(c) * sin(self.lat0) + y * sin(c) * cos(self.lat0) / p)
+    lon, lat = deg(lon), deg(lat)
+    return lon, lat
+end
+
+function AzimuthalEqualArea:model()
+    return {type="sphere", r=self.r}
+end
+
+local projs = {
+    AzimuthalEqualArea=AzimuthalEqualArea
+}
+
+-- Generic projection interface.
+
+local function Proj(name, origin, radius)
+    local proj = {}
+    proj.name = name
+    proj.lon0, proj.lat0 = unpack(origin or {0, 0})
+    proj.lon0, proj.lat0 = rad(proj.lon0), rad(proj.lat0)
+    proj.r = radius or 6378137
+    return setmetatable(proj, projs[name])
+end
+
+-- == Frames ==
+
+--[[
+ A frame stores information that specify the relation between the Earth's
+geometry (geodesy) and a particular map geometry (usually in 2D). Each map has a
+frame over which multiple layers of entities are drawn. There are two parameters
+that define a unique frame:
+* a projection;
+* a bounding box on the projected (plane) space.
+
+ The projection parameter must be fully specified, including the center of the
+projection, the orientation of the projection and the Earth model used (sphere
+or ellipsoid) along with its parameters.
+
+ The bounding box is specified as two points, determining minimal an maximal
+coordinates. The coordinate system for the bounding box is the projected one,
+but without scale, i.e. with meters as unit.
+]]
+
+local sep = ":"
+
+local function save_frame(fname, model, projection, bounding)
+    local frm = io.open(fname, "w")
+    frm:write("type", sep, model.type, "\n")
+    if model.type == "ellipsoid" then
+        frm:write("a", sep, model.a, "\n")
+        frm:write("b", sep, model.b, "\n")
+        frm:write("e", sep, model.e, "\n")
+        frm:write("f", sep, model.f, "\n")
+    elseif model.type == "sphere" then
+        frm:write("r", sep, model.r, "\n")
+    end
+    frm:write("proj", sep, projection.name, "\n")
+    frm:write("lon", sep, math.deg(projection.lon0), "\n")
+    frm:write("lat", sep, math.deg(projection.lat0), "\n")
+    frm:write("x0", sep, bounding.x0, "\n")
+    frm:write("y0", sep, bounding.y0, "\n")
+    frm:write("x1", sep, bounding.x1, "\n")
+    frm:write("y1", sep, bounding.y1, "\n")
+    frm:close()
+end
+
+local function load_frame(fname)
+    local frm = io.open(fname, "r")
+    local function get(field)
+        local line = frm:read()
+        local got = line:sub(1, #field)
+        assert(got == field, "expected field "..field.." but got "..got)
+        return line:sub(#field+#sep+1)
+    end
+    local model = {}
+    model.type = get "type"
+    if model.type == "ellipsoid" then
+        model.a = tonumber(get "a")
+        model.b = tonumber(get "b")
+        model.e = tonumber(get "e")
+        model.f = tonumber(get "f")
+    elseif model.type == "sphere" then
+        model.r = tonumber(get "r")
+    end
+    local projection = {}
+    projection.name = get "proj"
+    projection.lon = tonumber(get "lon")
+    projection.lat = tonumber(get "lat")
+    local bounding = {}
+    bounding.x0 = tonumber(get "x0")
+    bounding.y0 = tonumber(get "y0")
+    bounding.x1 = tonumber(get "x1")
+    bounding.y1 = tonumber(get "y1")
+    frm:close()
+    return model, projection, bounding
+end
+
+local Frame = {}
+Frame.__index = Frame
+
+function Frame:map(lon, lat)
+    local x, y = self.prj:map(lon, lat)
+    x = (x - self.bb.x0) * self.s
+    y = self.h - (y - self.bb.y0) * self.s
+    return x, y
+end
+
+function Frame:set_height(h)
+    local mw = self.bb.x1 - self.bb.x0
+    local mh = self.bb.y1 - self.bb.y0
+    self.h = h
+    self.s = h / mh
+    self.w = math.floor(mw * self.s + 0.5)
+end
+
+function Frame:mapped(points)
+    return coroutine.wrap(function()
+        for p in points do
+            local lat, lon = unpack(p)
+            local x, y = self:map(lat, lon)
+            coroutine.yield({x, y})
+        end
+    end)
+end
+
+local function new_frame(fname)
+    local self = setmetatable({}, Frame)
+    local m, p, b = load(fname)
+    self.bb = b
+    self.prj = proj.Proj(p.name, {p.lon, p.lat}, m.r)
+    return self
+end
+
+return {
+    distance=distance, bbox=bbox, centroid=centroid, Proj=Proj,
+    save_frame=save_frame, load_frame=load_frame, new_frame=new_frame
+}