#moongeodesy

中文简介见下方。

#Install

moon add ChrisTu117/moongeodesy

import {
  "ChrisTu117/moongeodesy/geodesic",
  "ChrisTu117/moongeodesy/utm",
  "ChrisTu117/moongeodesy/conic",
  "ChrisTu117/moongeodesy/magnetic",
  "ChrisTu117/moongeodesy/cn",
}

#Usage

let g = @geodesic.Geodesic::wgs84()

// Inverse problem: distance and bearings between two points
let r = g.inverse(39.9087, 116.3975, 31.2304, 121.4737) // Beijing → Shanghai
// r.s12 = 1066676.4 m, r.azi1 = 152.94°, r.azi2, r.m12, r.area12 ...

// Direct problem: go 1000 km due east from Shenzhen
let d = g.direct(22.5431, 114.0579, 90.0, 1000000.0) // d.lat2, d.lon2, d.azi2

// Geodesic line: waypoints at most 100 km apart, e.g. for drawing a route
let line = g.inverse_line(31.2304, 121.4737, 37.7749, -122.4194)
let pts = line.densify(100000.0)

// Polygon area (m², counter-clockwise positive) and perimeter (m)
let a = g.polygon_area([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)])

// Rhumb line: constant-azimuth distance, London → Tokyo (longer than the geodesic)
let rh = @geodesic.Rhumb::wgs84()
let h = rh.inverse(51.4778, -0.0015, 35.6762, 139.6503) // h.s12, h.azi12, h.area12

// Intersection of two geodesics (signed metres from each start)
let ix = @geodesic.Intersect::wgs84()
let cross = ix.closest(0.0, 0.0, 45.0, 45.0, 10.0, 135.0) // cross.x, cross.y, cross.c

// Local East-North-Up relative to an origin (ellipsoidal height in metres)
let lc = @geodesic.LocalCartesian::wgs84(39.9087, 116.3975)
let enu = lc.forward(39.9097, 116.3985, 50.0) // enu.x east, enu.y north, enu.z up
let back = lc.reverse(enu.x, enu.y, enu.z) // back.lat, back.lon, back.h

// Other ellipsoids: CGCS2000 / GRS80 built in, or any (a, f)
let cgcs = @geodesic.Geodesic::cgcs2000()
let sphere = @geodesic.Geodesic::new(6371000.0, 0.0)

// High flattening (|f| above about 1/50): elliptic-integral geodesic
let flat = @geodesic.GeodesicExact::new(6378137.0, 0.1)
let hx = flat.inverse(40.6, -73.8, 1.5, 0.0) // hx.s12

// Chinese map coordinate systems
let gcj = @cn.wgs84_to_gcj02(@cn.LatLon::new(39.9087, 116.3975)) // GPS → Gaode
let bd = @cn.gcj02_to_bd09(gcj) // Gaode → Baidu
let wgs = @cn.bd09_to_wgs84(bd) // exact inverse (round trip < 1e-9°)

// UTM / UPS / MGRS (functions raise @utm.UtmError on bad input)
let u = @utm.from_latlon(39.9087, 116.3975) // zone 50N, easting/northing, convergence, scale
let m = @utm.to_mgrs(39.9087, 116.3975) // "50SMK4850217797"
let p = @utm.from_mgrs("50S MK 48502 17797") // centre of the 1 m square

// Conic grids: CGCS2000 Albers (equal area) with parallels 25°N/47°N, 105°E
let albers = @conic.ConicGrid::albers(
  6378137.0,
  1.0 / 298.257222101,
  lat_origin=0.0,
  lon_origin=105.0,
  stdlat1=25.0,
  stdlat2=47.0,
)
let gp = albers.forward(39.9087, 116.3975) // gp.easting = 956117.8 m, gp.northing, gp.gamma, gp.k

// Magnetic declination (World Magnetic Model 2025) for a compass
let wmm = @magnetic.MagneticModel::wmm2025()
let t = @magnetic.decimal_year(2026, 11, 10)
let mag = wmm.field(39.9087, 116.3975, 0.0, t) // mag.decl = -7.56° (west), mag.incl, mag.h ...
let compass = @magnetic.magnetic_bearing(r.azi1, mag.decl) // true bearing → compass bearing

moon run src/examples/demo

moon run src/examples/cli -- geod inverse 39.9087 116.3975 31.2304 121.4737
moon run src/examples/cli -- geod direct 22.5431 114.0579 90 1000000
moon run src/examples/cli -- rhumb inverse 51.4778 -0.0015 35.6762 139.6503
moon run src/examples/cli -- intersect closest 0 0 45 45 10 135
moon run src/examples/cli -- enu forward 39.9087 116.3975 39.9097 116.3985 50
moon run src/examples/cli -- coord wgs2gcj 39.9087 116.3975
moon run src/examples/cli -- coord gcj2wgs 39.91010 116.40374
moon run src/examples/cli -- coord wgs2bd 39.9087 116.3975
moon run src/examples/cli -- coord bd2wgs 39.91645 116.41012
moon run src/examples/cli -- mag 39.9087 116.3975 2026-11-10
moon run src/examples/cli -- albers forward 25 47 0 105 39.9087 116.3975
moon run src/examples/cli -- lcc reverse 25 47 0 105 1000000 4000000

#Scenarios

#API reference

• geodesic — ellipsoids, inverse and direct, lines, polygons, rhumb lines, intersections, LocalCartesian ENU
• utm — transverse Mercator, UTM/UPS, MGRS
• conic — Lambert conformal conic, Albers equal-area, EPSG grids
• magnetic — World Magnetic Model 2025, bearings, grivation
• cn — WGS-84, GCJ-02, BD-09

moon doc --serve

#Accuracy and testing

#Notes

• GCJ-02 is the published reverse-engineered algorithm used by most open-source libraries (Krasovsky ellipsoid, a = 6378245 m). The transform is the identity outside the usual China bounding box, on the land of other countries inside it (Seoul, Hanoi, Naha, Penang, ...), and in the sea within about 9 km of their coasts, using Natural Earth 1:10m countries (China point of view) and minor islands. Islets missing from Natural Earth are treated as sea and offset; among GeoNames towns of 5000+ people there are 3 such cases (Kiên Hải, Vietnam; Lipayran and Litayan, Philippines). Within about 2 km of China's land border either answer is possible. Open-source implementations are reported to agree with the official, unpublished cipher to about 1 m; this library matches them to 1e-11° but cannot check the official one.
• The series geodesic (Geodesic) assumes |f| up to about 1/50 and the rhumb series about 1/150 (GeographicLib's guidance). GeodesicExact is the elliptic-integral solution for larger flattening, including prolate ellipsoids.
• Nearly antipodal points may have more than one shortest geodesic; then the distance is unique but azimuths/area depend on which one is returned (same behaviour as GeographicLib).
• A rhumb line that crosses a pole has no single longitude; lon2 and area12 are NaN in that case.
• Two closest geodesic intersections can be the same L1 distance apart. closest, next, and segment then return one of them. segment needs one shortest geodesic between each pair of endpoints.
• Karney (2011) shows that the transverse Mercator series is accurate to a few nanometres within ~3900 km of the central meridian; the tests here reach ~930 km.
• WMM2025 is valid from 2025.0 to 2030.0 and from 1 km below to 850 km above the ellipsoid; field raises outside that range. The tests check that the code evaluates the model correctly; the model's own declination uncertainty is about 0.3°–0.4° away from the magnetic poles.

#Benchmarks

#Port source and license

• geodesic (geodesics, GeodesicExact, polygons, rhumb lines, intersections, Geocentric and LocalCartesian), utm, and conic are ports of the algorithms in GeographicLib by Charles F. F. Karney (C++ and Python versions, v2.x), released under the MIT/X11 License, © Charles Karney. Rhumb lines follow GeographicLib's Rhumb class (order-6 series; Karney 2024, doi:10.1007/s11200-024-0709-z). GeodesicExact follows GeographicLib's elliptic-integral geodesic (Karney 2012–2023; area via the DST in GeographicLib 2.5). Intersections follow GeographicLib's Intersect class (Karney 2024, doi:10.1061/JSUED2.SUENG-1483; preprint arXiv:2308.00495). conic follows LambertConformalConic and AlbersEqualArea. Structure, series coefficients and edge-case handling follow GeographicLib; the MoonBit code is a rewrite, not a binding.
• magnetic is written from the equations in the WMM2025 technical report (Chulliat et al. 2025, doi:10.25923/prbc-s316), not ported from a library. The coefficients and test values are from NOAA NCEI and the British Geological Survey (doi:10.25921/aqfd-sd83); they are U.S. Government material, not subject to copyright protection (17 U.S.C. 403).
• GCJ-02/BD-09 use the widely published reverse-engineered formulas (as in coord-convert, MIT); the iterative exact inverses are our own.
• MGRS letter tables follow NGA's GeoTrans / GeographicLib MGRS.
• The China border data is Natural Earth 1:10m countries and minor islands (public domain).
• Test data: GeodTest (Karney, CC0); GeoNames towns (CC BY 4.0, in src/cn/towns_data_test.mbt); IOGP GIGS test dataset 2.1 (© IOGP, used under the GIGS Terms of Use, in src/conic/gigs_data_test.mbt only); the other vectors are generated by tools/ with geographiclib, PROJ (pyproj), mgrs, coord-convert, GeographicLib command-line tools, shapely, and mpmath.
• This library is MIT licensed (see LICENSE). The notices that third-party licences require, including GeographicLib's MIT notice, are in THIRD_PARTY_NOTICES.md.

#Code size

#History

#November 2026 contest period

• Fixes from the October acceptance review (#32): foreign land inside the China box is no longer offset, the exact GCJ-02 inverse is exact next to the box edges and borders, UPS works at ±90°.
• World Magnetic Model 2025, package magnetic, with CLI mag and scenario s5 (#33, #34).
• Lambert conformal conic and Albers equal-area projections, package conic, with CLI lcc / albers and scenario s6 (#35, #36).
• Foreign coasts and islands inside the box are no longer offset either: the boundary now uses Natural Earth minor islands and 9 km of coastal sea, and is tested on every GeoNames town of 5000+ people. Also: Albers poles round-trip exactly, and inverse conformal latitudes converge on prolate ellipsoids (0.2.1).
• Every README statement mapped to a test or check (docs/十一月验收清单.md): independent reference data for Geocentric/LocalCartesian, assertions tightened to the published bounds, scenarios that abort on wrong results and run on every backend in CI, and tests for API promises that had none.

#October 2026 contest period

#AI assistance

• Goal. A dependency-free MoonBit library for the geodesy that GIS, maps and logistics code needs in China: ellipsoidal distance/area, UTM/MGRS and conic grids, compass declination, and the GCJ-02/BD-09 offsets used by Gaode, Tencent and Baidu.
• Approach. An AI coding assistant translated GeographicLib's algorithms function by function, wrote the WMM2025 code from the NOAA report, and wrote tests; the author chose scope and package boundaries, reviewed the numerical code against the papers and sources, and decided on API shape and error handling. docs/十一月新增功能说明.md explains the November code in Chinese.
• How correctness is verified. Not by reading the AI's output, but by reference data that does not come from this code: Karney's official GeodTest data, PROJ, NGA GeoTrans, coord-convert, GeoNames, NOAA's WMM test values, IOGP GIGS, EPSG and Snyder worked examples, 40–50-digit mpmath evaluations of textbook formulas, and, for the GeographicLib ports, GeographicLib's own tools (table above). Every README statement has a test or check, listed in the acceptance checklist; CI runs everything on wasm, wasm-gc, js and native on Linux and macOS. The geodesic code needs a nearly correctly rounded hypot on very short lines; a named test checks two such lines against geographiclib.

#Related MoonBit libraries

• totto2727/geo — 2D planar geometry (port of georust/geo)
• cn-xjr/moongeokit — robust 2D geometry and spatial indexing
• LING71671/moon_strtree — STR-packed spatial indexes
• NBB2006/routeweave — route encoding, GeoJSON, Web Mercator tiles
• 2973181701/moongpx — GPX parsing
• liyun6666/moon-shapefile — Shapefile reader/writer

#Credits

• Geodesic and projection algorithms: C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 (2013), Transverse Mercator with an accuracy of a few nanometers, J. Geodesy 85, 475–485 (2011), and Geodesic intersections, J. Surveying Eng. 150(3), 04024005 (2024); code structure follows GeographicLib (MIT/X11 License, © Charles Karney).
• Conic projection formulas: IOGP Guidance Note 7-2 (EPSG methods 9801, 9802, 9822); J. P. Snyder, Map Projections — A Working Manual, USGS Professional Paper 1395 (1987).
• Geomagnetic model: A. Chulliat et al., The US/UK World Magnetic Model for 2025–2030: Technical Report, NOAA NCEI (2025).
• Test data: GeodTest.dat, Karney (2010), doi:10.5281/zenodo.32156, CC0; WMM2025 test values, NOAA NCEI; GIGS test dataset 2.1, © IOGP.

#License

#中文简介

• geodesic：移植 GeographicLib 的 Karney 测地线算法——椭球面上两点的最短距离与方位角（反算）、给定方位与距离求终点（正算）、测地线插点（用于地图画航线）、任意多边形的椭球面面积与周长、恒向线（rhumb line）、两条测地线的交点，以及相对原点的局部东北天直角坐标（LocalCartesian，经地心地固 ECEF）。级数解用于 |f| 约 1/50 以内（WGS-84 等）；|f| 更大时用 GeodesicExact（椭圆积分）。官方 GeodTest 测试集上反算距离误差 < 8 nm、正算位置误差 < 10 nm；恒向线与 RhumbSolve 2.5 对照，反算 < 4 nm、正算 < 11 nm；交点与 IntersectTool 2.5 对照 < 0.5 µm；ECEF 与局部 ENU 与 50 位精度的教科书公式及 PROJ 对照 < 20 nm。
• utm：横轴墨卡托投影（Krüger 6 阶级数）、极球面投影、UTM/UPS 坐标（含挪威、斯瓦尔巴特特殊分带）以及 MGRS 军事格网参考系的编码与解码。与 PROJ 结果相差 < 10 nm，MGRS 字符串与 NGA GeoTrans 完全一致。
• conic（11 月新增）：兰勃特等角圆锥投影和阿尔伯斯等积圆锥投影，单/双标准纬线，任意椭球，带 EPSG 假原点的格网（如 CGCS2000、25°N/47°N、中央经线 105°E 的全国 Albers 等积投影，用于面积统计）。与 50 位精度的 EPSG/Snyder 公式对照 < 25 nm，与 PROJ 对照 < 0.25 µm（两条标准纬线只差 1e-4° 时 PROJ 自身的公式丢失有效数字，差 6 mm），并通过 IOGP GIGS 官方测试（容差 3–5 cm，实际差 1–5 mm，与 PROJ 相同）。
• magnetic（11 月新增）：世界地磁模型 WMM2025，给出任意地点和日期的磁偏角、磁倾角、磁场分量及年变化率，罗盘方位与真方位、格网方位互换，以及磁极附近罗盘不可用的“禁用区/警告区”判断。与 NOAA 官方 100 个测试点、技术报告逐步算例对照，并与一个独立的 40 位精度实现对照 < 1e-10 nT。
• cn：国内地图坐标系转换——WGS-84（GPS）↔ GCJ-02（高德、腾讯，"火星坐标"）↔ BD-09（百度），提供迭代精确逆转换（往返误差 < 1e-9°）以及常用的单步近似公式。框内其他国家的陆地（首尔、河内、那霸、槟城等）及其沿海约 9 公里内不做偏移；用 GeoNames 全部 5000 人以上城镇检验，只有 3 个位于 Natural Earth 未收录小岛上的城镇例外。

let g = @geodesic.Geodesic::wgs84()
let r = g.inverse(39.9087, 116.3975, 31.2304, 121.4737) // 北京 → 上海：r.s12 ≈ 1066.68 km
let gcj = @cn.wgs84_to_gcj02(@cn.LatLon::new(39.9087, 116.3975)) // GPS → 高德
let bd = @cn.gcj02_to_bd09(gcj) // 高德 → 百度
let mgrs = @utm.to_mgrs(39.9087, 116.3975) // "50SMK4850217797"
let rh = @geodesic.Rhumb::wgs84()
let h = rh.inverse(51.4778, -0.0015, 35.6762, 139.6503) // 伦敦 → 东京，恒向线
let ix = @geodesic.Intersect::wgs84()
let p = ix.closest(0.0, 0.0, 45.0, 45.0, 10.0, 135.0) // 两条测地线的最近交点
let lc = @geodesic.LocalCartesian::wgs84(39.9087, 116.3975) // 以天安门为原点的局部 ENU
let enu = lc.forward(39.9097, 116.3985, 50.0) // enu.x 东, enu.y 北, enu.z 天
let wmm = @magnetic.MagneticModel::wmm2025()
let mag = wmm.field(39.9087, 116.3975, 0.0, 2026.86) // 北京磁偏角 mag.decl ≈ -7.56°（偏西）

moon run src/examples/cli -- geod inverse 39.9087 116.3975 31.2304 121.4737
moon run src/examples/cli -- rhumb inverse 51.4778 -0.0015 35.6762 139.6503
moon run src/examples/cli -- intersect closest 0 0 45 45 10 135
moon run src/examples/cli -- enu forward 39.9087 116.3975 39.9097 116.3985 50
moon run src/examples/cli -- coord wgs2gcj 39.9087 116.3975
moon run src/examples/cli -- mag 39.9087 116.3975 2026-11-10
moon run src/examples/cli -- albers forward 25 47 0 105 39.9087 116.3975