#moongeodesy

中文简介见下方。

#Install

moon add ChrisTu117/moongeodesy

import {
  "ChrisTu117/moongeodesy/geodesic",
  "ChrisTu117/moongeodesy/utm",
  "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

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

#Scenarios

#API reference

• geodesic — ellipsoids, inverse and direct, lines, polygons, rhumb lines, intersections, LocalCartesian ENU
• utm — transverse Mercator, UTM/UPS, MGRS
• 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); outside the rough mainland China bounding box the transform is the identity. Accuracy of the forward transform relative to the official one is ~1 m.
• The series geodesic (Geodesic) assumes |f| up to about 1/50. GeodesicExact is the elliptic-integral solution for larger flattening, including prolate ellipsoids. It matches GeodSolve -E to a few nanometres for |f| up to 0.5 in the tests above.
• 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. The rhumb series assume |f| up to about 1/150.
• 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.
• The transverse Mercator series is accurate to a few nanometres within ~3900 km of the central meridian.

#Benchmarks

#Port source and license

• geodesic (geodesics, GeodesicExact, polygons, rhumb lines, and intersections) and utm 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). Structure, series coefficients and edge-case handling follow GeographicLib; the MoonBit code is a rewrite, not a binding.
• 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.
• This library is MIT licensed (see LICENSE).

#Code size

#What is new in this contest period (October 2026)

#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 grids 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 and wrote tests; the author chose scope and package boundaries, reviewed the numerical code against Karney's papers and the GeographicLib sources, and decided on API shape and error handling.
• How correctness is verified. Not by reading the AI's output, but by independent references: Karney's official GeodTest data, Python geographiclib, PROJ 9, NGA GeoTrans and coord-convert (table above). The reference vectors are regenerated by scripts in tools/; CI runs every test on wasm, wasm-gc, js and native on Linux and macOS. Bugs found this way (e.g. the need for an accurate hypot on very short lines) were fixed and are covered by tests.

#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).
• Test data: GeodTest.dat, Karney (2010), doi:10.5281/zenodo.32156, CC0.

#License

#中文简介

• geodesic：移植 GeographicLib 的 Karney 测地线算法——椭球面上两点的最短距离与方位角（反算）、给定方位与距离求终点（正算）、测地线插点（用于地图画航线）、任意多边形的椭球面面积与周长、恒向线（等方位角的 rhumb line）、两条测地线的交点（最近点、下一点、线段），以及相对原点的局部东北天直角坐标（LocalCartesian，经地心地固 ECEF）。级数解用于 |f| 约 1/50 以内（WGS-84 等）。|f| 更大时用 GeodesicExact（椭圆积分，与 GeographicLib GeodSolve -E 对照，|f| 到 0.5 时距离最大误差 7.5 nm）。测地线精度为纳米级（官方 GeodTest 测试集上最大误差 < 10 nm）；恒向线与 GeographicLib RhumbSolve 2.5 对照，距离最大误差 3.73 nm；交点与 IntersectTool 2.5 对照，参数最大误差 0.45 µm；局部 ENU 与 GeographicLib / PROJ topocentric 对照，正算最大误差 < 1 nm。
• utm：横轴墨卡托投影（Krüger 6 阶级数）、极球面投影、UTM/UPS 坐标（含挪威、斯瓦尔巴特特殊分带）以及 MGRS 军事格网参考系的编码与解码。与 PROJ 结果相差 < 10 nm，MGRS 字符串与 NGA GeoTrans 完全一致。
• cn：国内地图坐标系转换——WGS-84（GPS）↔ GCJ-02（高德、腾讯，"火星坐标"）↔ BD-09（百度），提供迭代精确逆转换（往返误差 < 1e-9°）以及常用的单步近似公式。

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 天

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