fe-map: optional measurement along a sphere (shape conceded, radius free, haversine distance); land dots in the 3D view

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

femap/fit.py

@@ -94,3 +94,79 @@ def sphere_fit(X):
     u = (X - c) / d[:, None]; m = u.mean(0); m /= np.linalg.norm(m)
     return dict(radius=r, offset=float(np.median(np.abs(d - r) / r)),
                 wrap=float(np.degrees(np.arccos(np.clip(u @ m, -1, 1))).max()), centre=c.tolist())
+
+
+# ---- measured along a sphere: concedes the shape, leaves its size free ----
+def _arcs(P, R, a, b):
+    nrm = np.linalg.norm(P, axis=1); U = P / nrm[:, None]
+    u, v = U[a], U[b]; c = np.clip((u * v).sum(1), -1, 1); s = np.linalg.norm(np.cross(u, v), axis=1)
+    th = np.arctan2(s, c); s = np.maximum(s, 1e-9)
+    gu = -(v - c[:, None] * u) / (s * nrm[a])[:, None]; gv = -(u - c[:, None] * v) / (s * nrm[b])[:, None]
+    return R * th, th, gu, gv
+
+
+class SphereProblem(Problem):
+    """Cities held on a sphere of free radius; each distance measured along its surface (the haversine
+    arc, R x central angle). Stage 1 fits the radius with the cities; later stages hold it."""
+
+    def _pairs(self, P, R):
+        return _arcs(P, R, self.fa, self.fb), _arcs(P, R, self.la, self.lb)
+
+    def _apply(self, ga, gl, fa_, la_, R):
+        (dd, th, gu, gv), (ld, lth, lgu, lgv) = fa_, la_; n = self.n; G = np.zeros((n, 3))
+        for q in range(3):
+            G[:, q] = (np.bincount(self.fa, ga * R * gu[:, q], n) + np.bincount(self.fb, ga * R * gv[:, q], n)
+                       + np.bincount(self.la, gl * R * lgu[:, q], n) + np.bincount(self.lb, gl * R * lgv[:, q], n))
+        return G.ravel(), float((ga * th).sum() + (gl * lth).sum())
+
+    def arc_l2(self, wl=100.0):
+        def f(x):
+            P = x[:-1].reshape(self.n, 3); R = x[-1]; F_, L_ = self._pairs(P, R); dd, ld = F_[0], L_[0]
+            e = dd / self.fl - 1; ex = np.maximum(0, ld / self.lu - 1)
+            G, gR = self._apply(2 * e / self.fl, wl * 2 * ex / self.lu, F_, L_, R)
+            return (e ** 2).sum() + wl * (ex ** 2).sum(), np.r_[G, gR]
+        return f
+
+    def arc_sparse(self, R, lo, hi, w, wl):
+        a = self.fl * (1 - lo + MARGIN); b = self.fl * (1 + hi - MARGIN); ub = self.lu * (1 - MARGIN)
+        def f(x):
+            P = x.reshape(self.n, 3); F_, L_ = self._pairs(P, R); dd, ld = F_[0], L_[0]
+            eh = np.maximum(0, dd / b - 1); el = np.maximum(0, 1 - dd / a)
+            sh = np.sqrt(eh ** 2 + DELTA ** 2); sl = np.sqrt(el ** 2 + DELTA ** 2)
+            ex = np.maximum(0, ld / ub - 1); sx = np.sqrt(ex ** 2 + DELTA ** 2)
+            G, _ = self._apply(w * (eh / sh / b - el / sl / a), wl * ex / sx / ub, F_, L_, R)
+            return (w * (sh + sl - 2 * DELTA)).sum() + (wl * (sx - DELTA)).sum(), G
+        return f
+
+    def arc_broken(self, P, R, lo, hi):
+        F_, L_ = self._pairs(P, R); r = F_[0] / self.fl
+        return dict(long=int((r > 1 + hi).sum()), short=int((r < 1 - lo).sum()), limits=int((L_[0] > self.lu * 1.0005).sum()))
+
+    def arc_best_effort(self, X2, lo, hi, rounds=ROUNDS, maxiter=30000, eps=EPS, limit_w=LIMIT_W):
+        """From a flat layout laid on a spherical cap whose radius is a third of the layout's width."""
+        X2 = X2 - X2.mean(0); R0 = float(np.ptp(X2, axis=0).max()) / 3
+        x = np.r_[np.c_[X2, np.full(self.n, R0)].ravel(), R0]
+        x = minimize(self.arc_l2(), x, jac=True, method="L-BFGS-B", options={"maxiter": maxiter, "maxcor": 30}).x
+        P, R = x[:-1].reshape(self.n, 3), float(x[-1])
+        k = sum(self.arc_broken(P, R, lo, hi).values()); best = (k, P.copy())
+        w = np.ones(len(self.fl)); wl = np.full(len(self.lu), limit_w)
+        for _ in range(rounds):
+            P = minimize(self.arc_sparse(R, lo, hi, w, wl), P.ravel(), jac=True, method="L-BFGS-B",
+                         options={"maxiter": maxiter, "maxcor": 30}).x.reshape(self.n, 3)
+            k = sum(self.arc_broken(P, R, lo, hi).values())
+            if k < best[0]:
+                best = (k, P.copy())
+            F_, L_ = self._pairs(P, R); r = F_[0] / self.fl
+            e = np.maximum(0, r / (1 + hi) - 1) + np.maximum(0, 1 - r / (1 - lo)); l = np.maximum(0, L_[0] / self.lu - 1)
+            w = eps / (e + eps); wl = limit_w * eps / (l + eps)
+        P = best[1]; U = P / np.linalg.norm(P, axis=1)[:, None]
+        return R, U, self.arc_broken(P, R, lo, hi)
+
+
+def true_offsets(U, lat, lon, R):
+    """Distance of each solved city from its true place, after the orthogonal transformation that best
+    lines the true positions up with the solved ones (used to place coastlines and to check the fit)."""
+    p, l = np.radians(lat), np.radians(lon)
+    T = np.c_[np.cos(p) * np.cos(l), np.cos(p) * np.sin(l), np.sin(p)]
+    Us, _, Vt = np.linalg.svd(U.T @ T); Q = Us @ Vt
+    return R * np.arccos(np.clip(((T @ Q.T) * U).sum(1), -1, 1))