nx_geo_webmerc.nx source

1// nx_geo_webmerc.nx -- LIB: GEO-019b WEB-MERCATOR (EPSG:3857) slippy z/x/y tiles -- the web de-facto 2// standard (OSM/Google/Mapbox tile numbering). Hardening of GEO-019 (which is the geodetic grid). 3// 4// xtile = floor((lon+180) / 360 * 2^z) (linear, same as geodetic) 5// ytile = floor((1 - asinh(tan(lat)) / pi) / 2 * 2^z) (the Mercator projection) 6// with asinh(tan(lat)) = ln(tan(pi/4 + lat/2)) computed sovereignly via fx.nx CORDIC sin/cos + fx_log2 7// (ln x = log2 x * ln2). All Q16.16 fixed point -> deterministic; no float, no proj4, no API. 8// 9// This is last-mile interop (matching the worldwide web tile standard) per the ecosystem-only law -- 10// the projection itself is sovereign. Valid for |lat| < ~85.05 deg (the Mercator cutoff); beyond that 11// tan -> inf and we clamp. Convention: lat/lon microdegrees; out2 = [xtile, ytile]. license_tier: ORIGINAL 12import "fx.nx" 13import "nx_syscalls.nx" 14const MERC_MAGIC_180000000: i64 = 180000000 15const MERC_MAGIC_360000000: i64 = 360000000 16 17const MERC_LN2: i64 = 45426 // round(ln(2) * 65536) -- converts log2 -> ln in Q16.16 18 19func geo_webmerc_pow2(z: i64) -> i64 { var n: i64 = 1; var t: i64 = 0; while t < z { n = n * 2; t = t + 1 } return n } 20 21// (z, lat, lon) microdeg -> Web-Mercator tile [xtile, ytile] in out2, clamped to [0, 2^z-1]. 22func geo_webmerc_tile(z: i64, lat: i64, lon: i64, out2: *i64) -> i64 { 23 let n: i64 = geo_webmerc_pow2(z) 24 25 // x: linear in longitude. 26 var tx: i64 = (lon + MERC_MAGIC_180000000) * n / MERC_MAGIC_360000000 27 if tx < 0 { tx = 0 } 28 if tx >= n { tx = n - 1 } 29 30 // y: Mercator. phi (rad, Q16.16) = lat_microdeg * 2pi / 360e6 ; arg = pi/4 + phi/2. 31 let phi_q: i64 = lat * FX_TWO_PI / MERC_MAGIC_360000000 32 let arg: i64 = FX_PI / 4 + phi_q / 2 33 let s: i64 = fx_sin(arg) 34 let c: i64 = fx_cos(arg) 35 var ty: i64 = 0 36 if c == 0 { 37 ty = 0 // at the pole-ward cutoff; clamp below 38 } else { 39 var tanv: i64 = s * FX_ONE / c // tan in Q16.16 40 if tanv <= 0 { tanv = 1 } // log domain guard 41 let psi: i64 = fx_mul(fx_log2(tanv), MERC_LN2) // ln(tan(arg)) in Q16.16 radians 42 let psi_over_pi: i64 = psi * FX_ONE / FX_PI // Q16.16 43 let frac_y: i64 = (FX_ONE - psi_over_pi) / 2 // Q16.16 in [0,1] 44 ty = frac_y * n / FX_ONE 45 } 46 if ty < 0 { ty = 0 } 47 if ty >= n { ty = n - 1 } 48 49 out2[0] = tx 50 out2[1] = ty 51 return 0 52}