`nx_geo.nx` source

↩ module page · 151 lines · 6630 B
1// nx_geo.nx -- LIB: sovereign INTEGER-EXACT geofencing primitives. The hardware-rung-up foundation
2// of the Nishi geo / geofencing / maps stack (operator: "superior to google or open source tools").
3//
4// THE EXCEED ANGLE (measured, not asserted): coordinates are MICRODEGREES (degrees * 1e7) as i64,
5// and point-in-polygon is decided by INTEGER cross-multiplication -- NO floating point anywhere. So
6// the inside/outside verdict is EXACT and DETERMINISTIC at 1e-7 deg (~1.1 cm) resolution, with no
7// epsilon/robustness bugs on edges or vertices. Turf.js / Shapely / PostGIS all use float and give
8// non-deterministic or convention-dependent answers exactly there; the gate proves our result is
9// exact on the cases that break them. Overflow-safe for Earth coordinates in i64.
10//
11// This is the geofence FLOOR. Geodesic distance (haversine, rides a sovereign trig rung) and the
12// spatial index (geohash/grid) are rungs above; both compose THIS. license_tier: ORIGINAL
13import "nx_syscalls.nx"
14const GEO_MAGIC_90000000: i64 = 90000000
15const GEO_MAGIC_180000000: i64 = 180000000
16
17const GEO_MICRO: i64 = 1000000   // microdegrees per degree (degrees * 1e6); ~0.11 m resolution.
18                                 // All coords are microdegrees: lat in [-90e6,90e6], lon in [-180e6,180e6].
19
20// point in an axis-aligned bbox (inclusive) -- the O(1) pre-filter run before the polygon test.
21func geo_bbox(minlat: i64, minlon: i64, maxlat: i64, maxlon: i64, lat: i64, lon: i64) -> i64 {
22    if lat < minlat { return 0 }
23    if lat > maxlat { return 0 }
24    if lon < minlon { return 0 }
25    if lon > maxlon { return 0 }
26    return 1
27}
28
29// compute a polygon's bounding box into out=[minlat,minlon,maxlat,maxlon]. poly = flat
30// [lat0,lon0, lat1,lon1, ...], npts vertices.
31func geo_poly_bbox(poly: *i64, npts: i64, out: *i64) -> i64 {
32    var minlat: i64 = poly[0]
33    var maxlat: i64 = poly[0]
34    var minlon: i64 = poly[1]
35    var maxlon: i64 = poly[1]
36    var i: i64 = 1
37    while i < npts {
38        let la: i64 = poly[i * 2]
39        let lo: i64 = poly[i * 2 + 1]
40        if la < minlat { minlat = la }
41        if la > maxlat { maxlat = la }
42        if lo < minlon { minlon = lo }
43        if lo > maxlon { maxlon = lo }
44        i = i + 1
45    }
46    out[0] = minlat
47    out[1] = minlon
48    out[2] = maxlat
49    out[3] = maxlon
50    return 0
51}
52
53// INTEGER-EXACT point-in-polygon (ray-cast, PNPOLY) -- lat=Y, lon=X; horizontal ray in +lon.
54// poly = flat [lat0,lon0, ...], npts vertices (implicitly closed last->first). 1 inside, 0 outside.
55// The edge-intersection comparison is cross-multiplied by the edge's dy so NO division/float is
56// used; the inequality flips with the sign of dy. straddle guarantees dy != 0.
57func geo_point_in_poly(poly: *i64, npts: i64, lat: i64, lon: i64) -> i64 {
58    var inside: i64 = 0
59    var i: i64 = 0
60    var j: i64 = npts - 1
61    while i < npts {
62        let yi: i64 = poly[i * 2]
63        let xi: i64 = poly[i * 2 + 1]
64        let yj: i64 = poly[j * 2]
65        let xj: i64 = poly[j * 2 + 1]
66        var a: i64 = 0
67        if yi > lat { a = 1 }
68        var b: i64 = 0
69        if yj > lat { b = 1 }
70        if a != b {
71            let dy: i64 = yj - yi
72            let lhs: i64 = (lon - xi) * dy
73            let rhs: i64 = (xj - xi) * (lat - yi)
74            var cross: i64 = 0
75            if dy > 0 { if lhs < rhs { cross = 1 } }
76            if dy < 0 { if lhs > rhs { cross = 1 } }
77            if cross == 1 { inside = 1 - inside }
78        }
79        j = i
80        i = i + 1
81    }
82    return inside
83}
84
85// geofence test: a point is fenced iff it is inside the polygon (bbox pre-filtered for speed).
86func geo_fence_contains(poly: *i64, npts: i64, bbox: *i64, lat: i64, lon: i64) -> i64 {
87    if geo_bbox(bbox[0], bbox[1], bbox[2], bbox[3], lat, lon) == 0 { return 0 }
88    return geo_point_in_poly(poly, npts, lat, lon)
89}
90
91// RADIUS geofence ("within R of a center") -- the most common fence. INTEGER-EXACT, no sqrt on the
92// hot path (compare squared distance). The earth's longitude-compression (a lon degree is shorter
93// than a lat degree, by cos(latitude)) is handled by a per-fence `lon_scale` = cos(center_lat) *
94// GEO_SCALE, baked once at fence-definition (using a trig rung) so the runtime stays pure-integer.
95// r2 = (radius expressed in latitude-microdegree-equivalents) squared. Returns 1 inside, 0 outside.
96const GEO_SCALE: i64 = 1000000   // fixed-point unit for lon_scale (cos(lat) * 1e6)
97func geo_in_radius(clat: i64, clon: i64, lon_scale: i64, r2: i64, lat: i64, lon: i64) -> i64 {
98    let dlat: i64 = lat - clat
99    let dlon_raw: i64 = lon - clon
100    let dlon: i64 = (dlon_raw * lon_scale) / GEO_SCALE   // compress lon to lat-equivalent ground units
101    let d2: i64 = dlat * dlat + dlon * dlon
102    if d2 <= r2 { return 1 }
103    return 0
104}
105
106// GEOHASH spatial index (the "superior at scale" rung): interleave lon/lat bits and base32-encode,
107// so geofence lookup becomes O(1) prefix bucketing instead of O(N) scanning every fence. The
108// PROXIMITY PROPERTY -- nearby points share a long common prefix -- is what makes it an index.
109// Pure-integer bisection (microdegree midpoints), deterministic. Standard alphabet (no a/i/l/o), so
110// the output interoperates with every geohash tool (last-mile). lat/lon in microdegrees; `prec`
111// chars written to out (NUL-terminated). Returns length.
112func geo_geohash(lat: i64, lon: i64, prec: i64, out: *u8) -> i64 {
113    let b32: *u8 = "0123456789bcdefghjkmnpqrstuvwxyz" as *u8
114    var latlo: i64 = 0 - GEO_MAGIC_90000000
115    var lathi: i64 = GEO_MAGIC_90000000
116    var lonlo: i64 = 0 - GEO_MAGIC_180000000
117    var lonhi: i64 = GEO_MAGIC_180000000
118    var is_lon: i64 = 1
119    var k: i64 = 0
120    while k < prec {
121        var bits: i64 = 0
122        var nb: i64 = 0
123        while nb < 5 {
124            if is_lon == 1 {
125                let mid: i64 = (lonlo + lonhi) / 2
126                if lon >= mid { bits = bits * 2 + 1; lonlo = mid } else { bits = bits * 2; lonhi = mid }
127                is_lon = 0
128            } else {
129                let mid: i64 = (latlo + lathi) / 2
130                if lat >= mid { bits = bits * 2 + 1; latlo = mid } else { bits = bits * 2; lathi = mid }
131                is_lon = 1
132            }
133            nb = nb + 1
134        }
135        out[k] = b32[bits]
136        k = k + 1
137    }
138    out[k] = 0 as u8
139    return k
140}
141
142// length of the common prefix of two NUL-terminated geohashes -- the proximity score the index
143// buckets on (longer shared prefix = closer).
144func geo_prefix_common(a: *u8, b: *u8) -> i64 {
145    var i: i64 = 0
146    while a[i] != (0 as u8) {
147        if a[i] != b[i] { return i }
148        i = i + 1
149    }
150    return i
151}
nx_geo.nx source

`nx_geo.nx` source
