1package cells23import (4 "math"5 "strings"67 "github.com/golang/geo/r3"8 "github.com/golang/geo/s2"910 "github.com/FelineStateMachine/atlas/format/semconv"11)1213// The S2 side: the sphere divided by the six faces of a cube projected onto14// it, each face a quadtree of exactly-nested quadrilaterals addressed by15// Hilbert-curve tokens.16//17// Nothing spherical is written here. `github.com/golang/geo/s2` is the18// canonical Google library and the seam's s2js is a port of it, so the two19// lanes hold the same arithmetic by provenance rather than by agreement —20// which is the whole reason the server may answer a containment question at21// all instead of asking the seam for its answer. What this file adds is the22// two things the port has that the library does not: the bridge from world23// pixels to the sphere through the ground's declared flattening, and the level24// bookkeeping in which the root is a place.25//26// LEVELS. The port counts the root as 0 and the six faces as 1; S2 counts the27// faces as 0. The two bookkeepings meet in [S2Level] and nowhere else, exactly28// as they meet in `portLevel` in analysis/cellsystems/s2.ts.29//30// THE ROOT. The root is `""` in every system, and S2's token vocabulary has no31// spelling for it — `CellIDFromToken("")` answers a cell id that is neither32// the sphere nor an error. Every function below that could be handed the root33// answers for the whole ground instead of asking the library a question it34// cannot hear.3536// deepestS2Level is where the telescope stops: cells a few kilometres across37// on a Mars-sized body, tokens six characters long. It is the seam's floor38// (analysis/cellsystems/s2.ts) and is spelled here so that an address arriving39// from a request is held to the same floor the navigator would have produced.40const deepestS2Level = 104142// S2InputLength is how many characters the navigator's field accepts: a43// level-10 token is six, and the field takes no more (the port's44// `inputLength`). It is spelled beside the floor it follows from, because the45// two move together or the field starts accepting addresses the telescope46// would immediately clamp.47const S2InputLength = 64849// Mapping is a world's declared flattening, read once and inverted per point.50//51// A geodesic system addresses a sphere, and a packed position is a pixel on a52// raster. The mapping is the whole of what turns one into the other, and a53// world that declares none cannot be divided by S2 at all — which is the54// refusal [Applicable] makes on the session's behalf, before any of this is55// reached.56type Mapping struct{ window semconv.Equirect }5758// MappingOf reads the world's declared flattening out of its attributes, or59// refuses. It mirrors `geoMapping` in analysis/semconv/geometry.ts, including60// its one subtlety: the window is only read when `atlas.geometry.projection`61// says which flattening the numbers are for. A world carrying an equirect62// window under some other projection is declaring something this cannot read,63// and a plausible answer would be worse than none.64func MappingOf(attrs map[string]string) (Mapping, bool) {65 if attrs[semconv.KeyGeometryProjection] != semconv.ProjectionEquirect {66 return Mapping{}, false67 }68 window, declared, err := semconv.EquirectOf(attrs)69 if !declared || err != nil {70 return Mapping{}, false71 }72 return Mapping{window: window}, true73}7475// LatLng puts a point back on the sphere: world pixels in, true degrees out.76//77// The x, y are the application's own — y **negative-down**, the space a pin's78// position lives in — and the sign is flipped here, on the one line where the79// two conventions meet. [semconv.Equirect.Invert] reads y down from the top of80// the raster window, which is the convention the attributes are written in and81// the one the seam's mapping speaks (see the same flip at `leafOf` in82// analysis/cellsystems/s2.ts).83func (m Mapping) LatLng(x, y float64) (lat, lon float64) {84 return m.window.Invert(x, -y)85}8687// S2Held is the held-cell question for an S2 address: the point goes back on88// the sphere, becomes the leaf cell holding it, and the held cell is asked89// whether that leaf is one of its descendants.90//91// Two addresses hold nothing rather than something. An empty token is the92// root, and the root is the whole sphere — but the caller reaching this93// function has a cell held, and "no cell" is not a narrowing it should be94// asking about, so it is spelled out here as the whole ground rather than left95// to the library. An address that does not parse is not a place: nothing is96// standing in it, which is the conservative half of the disagreement — a map97// that draws nothing is wrong in a way a reader can see, where a server that98// narrows nothing while the seam narrows is wrong in a way nobody can.99func S2Held(mapping Mapping, token string) Held {100 if token == "" {101 return func(float64, float64) bool { return true }102 }103 cell, ok := s2Cell(token)104 if !ok {105 return func(float64, float64) bool { return false }106 }107 return func(x, y float64) bool {108 lat, lon := mapping.LatLng(x, y)109 return cell.Contains(s2.CellIDFromLatLng(s2.LatLngFromDegrees(lat, lon)))110 }111}112113// S2CellAt is the leaf-most address of a point at the telescope's floor: the114// answer the pin card carries and `locate` in the port. It is empty when the115// point is nowhere the mapping can name.116func S2CellAt(mapping Mapping, x, y float64) string {117 lat, lon := mapping.LatLng(x, y)118 leaf := s2.CellIDFromLatLng(s2.LatLngFromDegrees(lat, lon))119 if !leaf.IsValid() {120 return ""121 }122 return leaf.Parent(deepestS2Level).ToToken()123}124125// S2Level is the port level of an address: 0 for the root, 1 for a face, and126// so on down. A token that is not a place has no level, and says so with -1.127func S2Level(token string) int {128 if token == "" {129 return 0130 }131 cell, ok := s2Cell(token)132 if !ok {133 return -1134 }135 return cell.Level() + 1136}137138// S2Parent is the address one level out. The root has no parent and neither139// does an address nobody can read; both answer the root, which is where a140// reader telescoping out of an address they cannot be standing in belongs.141//142// This is the arithmetic that makes S2 the proof the contract holds nothing143// geohash-shaped: an S2 address does not refine by appending characters, so144// telescoping out of `47a1cb` is not dropping its last character — it is145// `47a1`, three characters shorter and spelled differently.146func S2Parent(token string) string {147 cell, ok := s2Cell(token)148 if !ok || cell.Level() == 0 {149 return ""150 }151 return cell.Parent(cell.Level() - 1).ToToken()152}153154// S2Children is the address's four children in their stable order, or the six155// faces when the address is the root. An address nobody can read has none.156func S2Children(token string) []string {157 if token == "" {158 out := make([]string, 0, 6)159 for face := range 6 {160 out = append(out, s2.CellIDFromFace(face).ToToken())161 }162 return out163 }164 cell, ok := s2Cell(token)165 if !ok {166 return nil167 }168 children := cell.Children()169 out := make([]string, 0, len(children))170 for _, child := range children {171 out = append(out, child.ToToken())172 }173 return out174}175176// S2ParseCell reads an address the way the navigator's field does: the root is177// a place, a token deeper than the telescope's floor clamps to the floor, and178// anything else is either a canonically spelled token or not a place at all.179// It is the port's `parseInput` (analysis/cellsystems/s2.ts) minus the180// keystroke normalization, which belongs to the field and not to the server.181func S2ParseCell(text string) (string, bool) {182 if text == "" {183 return "", true184 }185 cell, ok := s2Cell(text)186 if !ok {187 return "", false188 }189 if cell.Level() > deepestS2Level {190 return cell.Parent(deepestS2Level).ToToken(), true191 }192 return cell.ToToken(), true193}194195// S2Descend is the child of an address holding a point, or nothing at the196// telescope's floor. It is `descendTarget` in the port: the leaf under the197// point, taken back up to one level below the address given, which is the same198// child `contains` would have picked.199func S2Descend(mapping Mapping, token string, x, y float64) string {200 lat, lon := mapping.LatLng(x, y)201 leaf := s2.CellIDFromLatLng(s2.LatLngFromDegrees(lat, lon))202 if !leaf.IsValid() {203 return ""204 }205 target := 0206 if token != "" {207 cell, ok := s2Cell(token)208 if !ok {209 return ""210 }211 target = cell.Level() + 1212 }213 if target > deepestS2Level {214 return ""215 }216 return leaf.Parent(target).ToToken()217}218219// S2Center is a cell's centre as world pixels: the point on the sphere the220// library calls the middle of the cell, put back on the raster through the221// declared flattening. The root has no centre of its own and answers the222// middle of the ground it is standing for.223func S2Center(mapping Mapping, ground Extent, token string) (x, y float64) {224 cell, ok := s2Cell(token)225 if !ok {226 return (ground.MinX + ground.MaxX) / 2, (ground.MinY + ground.MaxY) / 2227 }228 return mapping.world(s2.CellFromCellID(cell).Center())229}230231// S2Extent is a cell's bounding rectangle in world pixels.232//233// It is a *drawn* rectangle rather than a spherical one, and the difference is234// the whole reason it is written out longhand. An S2 cell's edges are geodesics235// -- straight on no projection -- so the box a reader sees around one is the box236// around the flattened boundary, and it is that box the carry compares against a237// halved rectangle. A spherical area in steradians is not comparable with a238// rectangle in pixels at all, and a lat/lng bound would be a third number again:239// only the picture's own measure puts the two hierarchies on one scale.240//241// The walk is `ringOf` in analysis/cellsystems/s2.ts, step for step, including242// the sequential unwrap that keeps a loop crossing the antimeridian continuous243// and the closing point that joins the loop in the frame the walk ended in. A244// cell circling a pole accumulates a whole world of longitude, and both of those245// are what stop it from measuring the entire planet.246func S2Extent(mapping Mapping, ground Extent, token string) Extent {247 cell, ok := s2Cell(token)248 if !ok {249 return ground250 }251 points := s2Ring(mapping, ground, cell)252 out := Extent{MinX: math.Inf(1), MinY: math.Inf(1), MaxX: math.Inf(-1), MaxY: math.Inf(-1)}253 for _, at := range points {254 out.MinX, out.MaxX = math.Min(out.MinX, at[0]), math.Max(out.MaxX, at[0])255 out.MinY, out.MaxY = math.Min(out.MinY, at[1]), math.Max(out.MaxY, at[1])256 }257 return out258}259260// s2Ring walks the four geodesic edges in small great-circle steps and lands261// each one on the raster. The step count is the seam's: enough that a face262// reads as a curve and few enough that a cell at the floor is four corners.263func s2Ring(mapping Mapping, ground Extent, cell s2.CellID) [][2]float64 {264 face := s2.CellFromCellID(cell)265 steps := int(math.Max(2, math.Min(16, math.Pow(2, float64(4-cell.Level()))*2)))266 width := ground.MaxX - ground.MinX267 points := make([][2]float64, 0, 4*steps+1)268 for edge := range 4 {269 from := face.Vertex(edge).Vector270 to := face.Vertex((edge + 1) % 4).Vector271 for step := range steps {272 part := float64(step) / float64(steps)273 at := r3.Vector{274 X: from.X + (to.X-from.X)*part,275 Y: from.Y + (to.Y-from.Y)*part,276 Z: from.Z + (to.Z-from.Z)*part,277 }278 norm := math.Sqrt(at.X*at.X + at.Y*at.Y + at.Z*at.Z)279 if norm == 0 {280 norm = 1281 }282 x, y := mapping.world(s2.Point{Vector: r3.Vector{X: at.X / norm, Y: at.Y / norm, Z: at.Z / norm}})283 points = append(points, [2]float64{x, y})284 }285 }286 // Sequential unwrap: a jump of more than half the world is the seam, not287 // the ground.288 for at := 1; at < len(points); at++ {289 delta := points[at][0] - points[at-1][0]290 if delta > width/2 {291 points[at][0] -= width292 } else if delta < -width/2 {293 points[at][0] += width294 }295 }296 closing := points[0]297 delta := closing[0] - points[len(points)-1][0]298 if delta > width/2 {299 closing[0] -= width300 } else if delta < -width/2 {301 closing[0] += width302 }303 return append(points, closing)304}305306// world is a point on the sphere as world pixels, y negative-down: the reverse307// of [Mapping.LatLng] and the same one flip of the sign.308func (m Mapping) world(point s2.Point) (x, y float64) {309 at := s2.LatLngFromPoint(point)310 x, y = m.window.Apply(at.Lat.Degrees(), at.Lng.Degrees())311 return x, -y312}313314// s2Cell reads a token, refusing what is not a cell.315//316// The guard is not a formality. A token the library cannot read answers the317// zero cell id, whose descendant range is the whole of the space — so an318// unguarded `Contains` on it says yes to every point on the sphere, and a319// misspelled address would quietly stop narrowing anything.320func s2Cell(token string) (s2.CellID, bool) {321 cell := s2.CellIDFromToken(strings.TrimSpace(token))322 if cell == 0 || !cell.IsValid() {323 return 0, false324 }325 return cell, true326}