`nx_geo_join.nx` source

1// nx_geo_join.nx -- LIB: GEO-013 SPATIAL JOIN (point-in-which-polygon), INTEGER-EXACT. The PostGIS
2// ST_Contains-join staple: assign each of N query points to the polygon (from a pool of M) that
3// contains it. Composes the GEO-001 exact ray-cast predicate, here as a pooled, base-offset variant
4// `geo_pip_off` so many polygons share one flat vertex buffer (pool + per-polygon offset/count).
5//
6// EXCEED ANGLE (measured, not asserted): the containment test is the same integer cross-multiplied
7// ray-cast as GEO-001 -- ZERO floating point -> the point->polygon assignment is EXACT and
8// DETERMINISTIC on edges/vertices, where float spatial joins (Turf/Shapely/PostGIS) can drop or
9// double-count boundary points. Overflow-safe for Earth coords in i64.
10//
11// (geo_pip_off restates GEO-001's ray-cast for the pooled layout; if a 3rd caller needs it, extract a
12// shared base-offset primitive into nx_geo per the DRY rule.) license_tier: ORIGINAL
13import "nx_syscalls.nx"
14
15// point-in-polygon for a polygon stored at `pool[(base+k)*2 .. ]`, npts vertices. 1 inside, 0 outside.
16// Same integer ray-cast as GEO-001 geo_point_in_poly (cross-multiplied by edge dy; no division/float).
17func geo_pip_off(pool: *i64, base: i64, npts: i64, lat: i64, lon: i64) -> i64 {
18    var inside: i64 = 0
19    var i: i64 = 0
20    var j: i64 = npts - 1
21    while i < npts {
22        let yi: i64 = pool[(base + i) * 2]
23        let xi: i64 = pool[(base + i) * 2 + 1]
24        let yj: i64 = pool[(base + j) * 2]
25        let xj: i64 = pool[(base + j) * 2 + 1]
26        var a: i64 = 0
27        if yi > lat { a = 1 }
28        var b: i64 = 0
29        if yj > lat { b = 1 }
30        if a != b {
31            let dy: i64 = yj - yi
32            let lhs: i64 = (lon - xi) * dy
33            let rhs: i64 = (xj - xi) * (lat - yi)
34            var cross: i64 = 0
35            if dy > 0 { if lhs < rhs { cross = 1 } }
36            if dy < 0 { if lhs > rhs { cross = 1 } }
37            if cross == 1 { inside = 1 - inside }
38        }
39        j = i
40        i = i + 1
41    }
42    return inside
43}
44
45// index of the first polygon (of npolys) containing (lat,lon); -1 if none. off[p]/cnt[p] = polygon p's
46// base vertex offset and vertex count within the shared pool.
47func geo_locate(pool: *i64, off: *i64, cnt: *i64, npolys: i64, lat: i64, lon: i64) -> i64 {
48    var p: i64 = 0
49    while p < npolys {
50        if geo_pip_off(pool, off[p], cnt[p], lat, lon) == 1 { return p }
51        p = p + 1
52    }
53    return 0 - 1
54}
55
56// spatial join: for each of npts query points write the containing polygon index (-1 if none) to out.
57func geo_spatial_join(pts: *i64, npts: i64, pool: *i64, off: *i64, cnt: *i64, npolys: i64, out: *i64) -> i64 {
58    var i: i64 = 0
59    while i < npts {
60        out[i] = geo_locate(pool, off, cnt, npolys, pts[i * 2], pts[i * 2 + 1])
61        i = i + 1
62    }
63    return npts
64}
nx_geo_join.nx source

`nx_geo_join.nx` source
