nx_geo_tiles.nx source

1// nx_geo_tiles.nx -- LIB: GEO-019 VECTOR/RASTER TILE addressing (sovereign slippy z/x/y scheme). 2// 3// Uses the GEODETIC / plate-carree tiling (the OGC WMTS WGS84 grid): the world lon[-180,180] x 4// lat[-90,90] is split into 2^z x 2^z tiles, LINEARLY in lon/lat -- so the whole scheme is INTEGER: 5// tile_x = floor((lon+180e6) * 2^z / 360e6), tile_y = floor((90e6-lat) * 2^z / 180e6) (y south+) 6// Honest: this is the geodetic grid, NOT Web-Mercator (Mercator needs log/tan -- a later rung). 7// 8// THE EXCEED ANGLE (measured): tile boundaries are exact integer microdegree edges, so tiles tessellate 9// the plane with NO seam gaps or overlaps -- the round-trip "a point's tile's bbox contains the point" 10// holds exactly (for z<=8, where 2^z divides 360e6 and 180e6; larger z floors the edges). Float tile 11// math (Mapbox/Google JS) can produce sub-pixel seam gaps/overlaps at boundaries. license_tier: ORIGINAL 12import "nx_syscalls.nx" 13 14const TILE_LON_SPAN: i64 = 360000000 15const TILE_LAT_SPAN: i64 = 180000000 16const TILE_LON_MIN: i64 = 0 - 180000000 17const TILE_LAT_MAX: i64 = 90000000 18 19func geo_tile_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) -> tile (tx,ty) written to out2, clamped to [0, 2^z-1]. 22func geo_point_tile(z: i64, lat: i64, lon: i64, out2: *i64) -> i64 { 23 let n: i64 = geo_tile_pow2(z) 24 var tx: i64 = (lon - TILE_LON_MIN) * n / TILE_LON_SPAN 25 var ty: i64 = (TILE_LAT_MAX - lat) * n / TILE_LAT_SPAN 26 if tx < 0 { tx = 0 } 27 if tx >= n { tx = n - 1 } 28 if ty < 0 { ty = 0 } 29 if ty >= n { ty = n - 1 } 30 out2[0] = tx 31 out2[1] = ty 32 return 0 33} 34 35// (z, tx, ty) -> bbox [minlat, minlon, maxlat, maxlon] written to out4 (microdegrees). 36func geo_tile_bbox(z: i64, tx: i64, ty: i64, out4: *i64) -> i64 { 37 let n: i64 = geo_tile_pow2(z) 38 out4[0] = TILE_LAT_MAX - (ty + 1) * TILE_LAT_SPAN / n // minlat (south edge) 39 out4[1] = TILE_LON_MIN + tx * TILE_LON_SPAN / n // minlon (west edge) 40 out4[2] = TILE_LAT_MAX - ty * TILE_LAT_SPAN / n // maxlat (north edge) 41 out4[3] = TILE_LON_MIN + (tx + 1) * TILE_LON_SPAN / n // maxlon (east edge) 42 return 0 43} 44 45// 1 iff (lat,lon) lies within the bbox out4 (inclusive). The round-trip containment predicate. 46func geo_bbox_contains(bbox: *i64, lat: i64, lon: i64) -> i64 { 47 if lat < bbox[0] { return 0 } 48 if lat > bbox[2] { return 0 } 49 if lon < bbox[1] { return 0 } 50 if lon > bbox[3] { return 0 } 51 return 1 52}