National hydrography enrichment
- commit
- 81d35937855b1c811b20c19f1bc21369ea99952b
- parent
- 68b174f5ad3b39f611df4c40e50e96c8fd1695cf
- parent
- 9056ca029a102273c57d44a64412c09033fdf5e7
- author
- FelineStateMachine <dami.etoile@gmail.com>
- kind
- feat risk score 2.1 of 3, Jev's estimate
| M | docs/ | 4 | + |
| M | internal/ | 48 | +++++- |
| M | internal/ | 85 | ++++++++++- |
| A | internal/ | 203 | +++++++++++++++++++++++++ |
| A | internal/ | 144 | ++++++++++++++++++ |
| A | internal/ | 240 | ++++++++++++++++++++++++++++++ |
| A | internal/ | 106 | +++++++++++++ |
| A | internal/ | 108 | ++++++++++++++ |
| M | internal/ | 1 | + |
| M | internal/ | 20 | +++ |
| M | internal/ | 2 | +- |
| M | tools/ | 26 | ++- |
| A | tools/ | 193 | ++++++++++++++++++++++++ |
| A | tools/ | 120 | +++++++++++++++ |
@@ -13,6 +13,10 @@ in a repository of their own, so every employee works off the same enriched
1313 1414 1. Crawls the city into the `fmg-archive` committed on the state branch.1515 Captures are content-addressed: an unchanged city registers nothing.16+ Beside the city's own hub, every crawl queries the USGS national17+ hydrography services on `hydro.nationalmap.gov` for the city's window —18+ watersheds, subwatersheds, named streams, waterbodies — so a runner19+ with an egress allowlist needs that host beside the hub's.1620 2. Composes the bundle. Its filename embeds the capture day and a content1721 stamp over everything inside.1822 3. Gates: if that exact filename is already a release asset, the run ends —@@ -88,6 +88,17 @@ type Dataset struct {
8888 // than a lane -- and absent means every feature draws alike.8989 Role string9090 Emphasis func(Fields) float6491+92+ // National enrichment. A dataset with Server set is fetched from a93+ // plain ArcGIS MapServer query endpoint by the city's padded window94+ // rather than from the city's hub, Where narrowing the rows, and its95+ // geometry is cut to the window at capture time. The first kept field96+ // names the layer's row identity and orders the pages; IDOf derives97+ // the capture's stable feature identifier from the kept fields when98+ // the layer's own object ids churn between upstream refreshes.99+ Server string100+ Where string101+ IDOf func(Fields) (int64, bool)91102 }92103 93104 // City is one curated hub. The table is the whole authority: an uncurated@@ -103,6 +114,13 @@ type City struct {
103114 Datasets []Dataset104115 }105116 117+// AllDatasets is the city's own curation followed by the national118+// enrichment every US city receives by bounding box, in that order, so a119+// city's zones keep their places and the nation's sort after them.120+func (c City) AllDatasets() []Dataset {121+ return append(append(make([]Dataset, 0, len(c.Datasets)+len(National)), c.Datasets...), National...)122+}123+106124 // Cities is every city the source knows how to capture. A city that should107125 // stay out of the repository -- your own, say -- can be registered from a108126 // cities_local.go beside this file: that name is ignored by git, and an@@ -529,7 +547,7 @@ func Translate(doc []byte) ([]byte, error) {
529547 if err != nil {530548 return nil, err531549 }532- regions, err := buildRegions(&capture, byOrder, ids)550+ regions, err := buildRegions(&capture, byOrder, ids, buildHydroIndex(&capture))533551 if err != nil {534552 return nil, err535553 }@@ -590,8 +608,9 @@ func curatedOrder(city City, capture *Capture) ([]pairing, error) {
590608 captured[dataset.Slug] = dataset591609 }592610 var out []pairing593- for at := range city.Datasets {594- curated := &city.Datasets[at]611+ curatedAll := city.AllDatasets()612+ for at := range curatedAll {613+ curated := &curatedAll[at]595614 if data, taken := captured[curated.Slug]; taken {596615 out = append(out, pairing{curated, data})597616 delete(captured, curated.Slug)@@ -739,8 +758,10 @@ const zoneLimit = 256
739758 // buildRegions folds the polygon and line datasets into zones: every740759 // feature whose ZoneOf names a bucket lands its ground in that bucket's741760 // region, and line features widen into ribbon polygons first, because a742-// zone is ground and a line has none.743-func buildRegions(capture *Capture, pairs []pairing, ids *mgdoc.IDSpace) ([]mgdoc.Region, error) {761+// zone is ground and a line has none. Zones of the city's own datasets --762+// never the national ones -- learn their subwatershed from the hydro763+// index when their features agree on one.764+func buildRegions(capture *Capture, pairs []pairing, ids *mgdoc.IDSpace, hydro *hydroIndex) ([]mgdoc.Region, error) {744765 regions := []mgdoc.Region{}745766 for _, pair := range pairs {746767 curated, data := pair.curated, pair.data@@ -748,12 +769,21 @@ func buildRegions(capture *Capture, pairs []pairing, ids *mgdoc.IDSpace) ([]mgdo
748769 continue749770 }750771 buckets := make(map[string]*mgdoc.Region)772+ claims := make(map[string]*hydroClaims)751773 var order []string752774 for _, feature := range data.Features {753775 key := curated.ZoneOf(feature.Fields)754776 if key.Key == "" {755777 continue756778 }779+ if hydro != nil && curated.Server == "" {780+ claim, tracked := claims[key.Key]781+ if !tracked {782+ claim = &hydroClaims{}783+ claims[key.Key] = claim784+ }785+ claim.observe(hydro, feature.Geometry)786+ }757787 zone, made := buckets[key.Key]758788 if !made {759789 if len(buckets) >= zoneLimit {@@ -788,9 +818,13 @@ func buildRegions(capture *Capture, pairs []pairing, ids *mgdoc.IDSpace) ([]mgdo
788818 }789819 sort.Strings(order)790820 for _, key := range order {791- if len(buckets[key].Features) > 0 {792- regions = append(regions, *buckets[key])821+ if len(buckets[key].Features) == 0 {822+ continue823+ }824+ if claim := claims[key]; claim != nil {825+ claim.apply(hydro, buckets[key])793826 }827+ regions = append(regions, *buckets[key])794828 }795829 }796830 return regions, nil@@ -104,6 +104,9 @@ func fixtureCapture() Capture {
104104 {ID: 1, Fields: Fields{"ZONE": "PF"}, Geometry: square(-121.31, 44.05, 0.004)},105105 {ID: 2, Fields: Fields{"ZONE": "RS"}, Geometry: square(-121.30, 44.06, 0.002)},106106 {ID: 3, Fields: Fields{"ZONE": "PF"}, Geometry: square(-121.28, 44.08, 0.003)},107+ // A second RS parcel in the other subwatershed: the RS zone108+ // straddles, and must say nothing about membership.109+ {ID: 4, Fields: Fields{"ZONE": "RS"}, Geometry: square(-121.36, 44.135, 0.002)},107110 }},108111 {Slug: "annexations", Features: []Feature{109112 {ID: 3, Fields: Fields{"DESCRIPT": "Rockridge Park", "EFF_DATE": "1556866800000", "ORDIN_NO": "NS-2327"}, Geometry: square(-121.29, 44.09, 0.002)},@@ -112,12 +115,31 @@ func fixtureCapture() Capture {
112115 {ID: 1, Fields: Fields{"TYPE": "Significant Identified Wetlands", "MAP_CODE": "R9"}, Geometry: square(-121.32, 44.04, 0.003)},113116 }},114117 {Slug: "trails", Features: []Feature{115- {ID: 6, Fields: Fields{"Trail_Name": "Riley Ranch", "Park": "Riley Ranch Nature Reserve", "Status": "Existing"}, Geometry: line([2]float64{-121.35, 44.11}, [2]float64{-121.345, 44.115})},118+ {ID: 6, Fields: Fields{"Trail_Name": "Riley Ranch", "Park": "Riley Ranch Nature Reserve", "Status": "Existing"}, Geometry: line([2]float64{-121.35, 44.14}, [2]float64{-121.345, 44.145})},116119 {ID: 7, Fields: Fields{"Trail_Name": "Planned Spur", "Status": "Proposed"}, Geometry: line([2]float64{-121.33, 44.02}, [2]float64{-121.329, 44.021})},117120 }},118121 {Slug: "mpo-boundary", Features: []Feature{119122 {ID: 21, Fields: Fields{"LABEL": "MPO BOUNDARY"}, Geometry: square(-121.40, 43.97, 0.17)},120123 }},124+ // The national enrichment, city-sized: two subwatersheds that125+ // between them hold most of the city fixtures -- but not the MPO126+ // boundary's corners, so that zone's ground lands nowhere and it127+ // claims nothing.128+ {Slug: "watersheds", Features: []Feature{129+ {ID: 1707030101, Fields: Fields{"huc10": "1707030101", "name": "Tumalo Creek"}, Geometry: square(-121.40, 44.00, 0.15)},130+ }},131+ {Slug: "subwatersheds", Features: []Feature{132+ {ID: 170703010101, Fields: Fields{"huc12": "170703010101", "name": "Tumalo Creek"}, Geometry: square(-121.37, 44.03, 0.10)},133+ {ID: 170703010102, Fields: Fields{"huc12": "170703010102", "name": "Bull Creek"}, Geometry: square(-121.37, 44.131, 0.03)},134+ }},135+ {Slug: "streams", Features: []Feature{136+ {ID: 51, Fields: Fields{"gnis_name": "Tumalo Creek"}, Geometry: line([2]float64{-121.34, 44.05}, [2]float64{-121.30, 44.07})},137+ }},138+ {Slug: "waterbodies", Features: []Feature{139+ {ID: 61, Fields: Fields{"GNIS_NAME": "Mirror Pond"}, Geometry: square(-121.316, 44.057, 0.002)},140+ // An unnamed pond draws into the basemap but earns no zone.141+ {ID: 62, Geometry: square(-121.29, 44.05, 0.001)},142+ }},121143 },122144 }123145 }@@ -204,13 +226,14 @@ func TestTranslateShapesTheDocument(t *testing.T) {
204226 205227 // Zones: the boundary as one, zoning bucketed by code, annexations by206228 // decade, wetlands by type, and only the existing named trail -- in207- // curated dataset order.229+ // curated dataset order, with the national enrichment after them all.208230 titles := make([]string, 0, len(out.Regions))209231 for _, region := range out.Regions {210232 titles = append(titles, region.Title)211233 }212234 want := []string{"MPO Boundary", "PF", "RS", "Annexed 2010–2019",213- "Significant Identified Wetlands", "Riley Ranch"}235+ "Significant Identified Wetlands", "Riley Ranch",236+ "Tumalo Creek", "Tumalo Creek", "Bull Creek", "Tumalo Creek", "Mirror Pond"}214237 if len(titles) != len(want) {215238 t.Fatalf("zones are %v", titles)216239 }@@ -219,6 +242,48 @@ func TestTranslateShapesTheDocument(t *testing.T) {
219242 t.Fatalf("zones are %v, want %v", titles, want)220243 }221244 }245+246+ // The national zones say which grain of the nation they are.247+ if got := out.Regions[6].Subtitle; got != "Watershed · HUC 1707030101" {248+ t.Fatalf("watershed subtitle is %q", got)249+ }250+ if got := out.Regions[7].Subtitle; got != "Subwatershed · HUC 170703010101" {251+ t.Fatalf("subwatershed subtitle is %q", got)252+ }253+ if got := out.Regions[9].Subtitle; got != "Stream" {254+ t.Fatalf("stream subtitle is %q", got)255+ }256+ if got := out.Regions[10].Subtitle; got != "Waterbody" {257+ t.Fatalf("waterbody subtitle is %q", got)258+ }259+260+ // Membership: the PF zone's ground lies wholly in one subwatershed, so261+ // it says so, once, in both spellings.262+ pf := out.Regions[1]263+ if pf.Description != "Lies in the Tumalo Creek subwatershed (HUC 170703010101)." {264+ t.Fatalf("PF card reads %q", pf.Description)265+ }266+ if err := semconv.Validate(semconv.EntityZone, pf.Attrs); err != nil {267+ t.Fatalf("PF attrs: %v", err)268+ }269+ if pf.Attrs[semconv.KeyHydroHUC12] != "170703010101" {270+ t.Fatalf("PF claims %q", pf.Attrs[semconv.KeyHydroHUC12])271+ }272+ // The RS zone straddles the two subwatersheds and says nothing; the MPO273+ // boundary's corners land in no unit and it says nothing either; the274+ // national zones never describe themselves.275+ for _, at := range []int{2, 0, 7} {276+ zone := out.Regions[at]277+ if zone.Description != "" || zone.Attrs[semconv.KeyHydroHUC12] != "" {278+ t.Fatalf("zone %q claims membership: %q %v", zone.Title, zone.Description, zone.Attrs)279+ }280+ }281+ // A line zone in one subwatershed claims it alongside its stroke.282+ riley := out.Regions[5]283+ if riley.Attrs[semconv.KeyHydroHUC12] != "170703010102" ||284+ riley.Attrs[semconv.KeyStrokeWidthPx] != "12" {285+ t.Fatalf("trail attrs are %v", riley.Attrs)286+ }222287 if out.Regions[1].Subtitle != "Zoning" {223288 t.Fatalf("PF subtitle is %q", out.Regions[1].Subtitle)224289 }@@ -246,7 +311,7 @@ func TestTranslateShapesTheDocument(t *testing.T) {
246311 if len(lineCoords) != 1 || len(lineCoords[0]) != 2 {247312 t.Fatalf("the two-point trail should stay one two-point line, got %v", lineCoords)248313 }249- if out.Regions[1].Attrs != nil {314+ if _, stroked := out.Regions[1].Attrs[semconv.KeyStrokeWidthPx]; stroked {250315 t.Fatalf("a polygon zone declares no stroke, got %v", out.Regions[1].Attrs)251316 }252317 @@ -426,6 +491,9 @@ func TestRound7(t *testing.T) {
426491 427492 // The curated table itself must be internally sound: slugs unique and428493 // well-formed, point datasets fully spelled, zone functions only on ground.494+// Every city walks with the national enrichment appended, so a private495+// city registered from cities_local.go is held to the same rules,496+// collisions with the national slugs included.429497 func TestCuratedTableIsSound(t *testing.T) {430498 for slug, city := range Cities {431499 if slug != city.Slug {@@ -435,13 +503,16 @@ func TestCuratedTableIsSound(t *testing.T) {
435503 t.Fatalf("city %q declares no ground or no pyramid", slug)436504 }437505 seen := map[string]bool{}438- for _, dataset := range city.Datasets {506+ for _, dataset := range city.AllDatasets() {439507 if dataset.Slug == "" || seen[dataset.Slug] {440508 t.Fatalf("dataset slug %q is empty or doubled", dataset.Slug)441509 }442510 seen[dataset.Slug] = true443- if dataset.ItemID == "" || dataset.Title == "" {444- t.Fatalf("dataset %q is missing its identity", dataset.Slug)511+ if (dataset.ItemID == "") == (dataset.Server == "") || dataset.Title == "" {512+ t.Fatalf("dataset %q wants exactly one identity, hub or national", dataset.Slug)513+ }514+ if dataset.Server != "" && len(dataset.Keep) == 0 {515+ t.Fatalf("dataset %q pages a national layer without a row identity", dataset.Slug)445516 }446517 isPoint := dataset.Geometry == "point"447518 if isPoint != (dataset.Group != "") {@@ -0,0 +1,203 @@
1+package arcgismap2+3+// Clipping for national datasets. A watershed intersecting a city's window4+// may extend a hundred windows past it, and the synthetic projection is5+// linear in the window's own terms: a vertex far outside would land on an6+// absurd longitude and wrap the viewer's world. City data never needs this7+// -- a hub publishes ground near its own city -- so only features fetched by8+// window are cut down to it, at capture time, where the trimmed shape rides9+// the content addressing like any other fact about the day.10+//11+// The window rectangle is cut in plain degrees. Its edges are lines of12+// constant longitude and latitude, and both the Mercator forward and the13+// synthetic spelling are monotone per axis, so the degree rectangle and the14+// projected square bound exactly the same ground.15+16+// ClipRings cuts MultiPolygon-nested rings to the window, ring by ring with17+// Sutherland–Hodgman. A polygon whose outer ring is cut to nothing is18+// dropped whole; a hole cut to nothing drops alone. Introduced vertices are19+// rounded like every captured coordinate, and output rings close on their20+// first position.21+func ClipRings(w Window, polygons [][][][]float64) [][][][]float64 {22+ var out [][][][]float6423+ for _, polygon := range polygons {24+ var kept [][][]float6425+ for at, ring := range polygon {26+ clipped := clipRing(w, ring)27+ if clipped == nil {28+ if at == 0 {29+ kept = nil30+ break31+ }32+ continue33+ }34+ kept = append(kept, clipped)35+ }36+ if len(kept) > 0 {37+ out = append(out, kept)38+ }39+ }40+ return out41+}42+43+// ClipLines cuts MultiLineString-nested lines to the window. A line that44+// leaves and returns splits into one part per crossing, and parts too short45+// to draw are dropped.46+func ClipLines(w Window, lines [][][]float64) [][][]float64 {47+ var out [][][]float6448+ var part [][]float6449+ flush := func() {50+ if len(part) >= 2 {51+ out = append(out, part)52+ }53+ part = nil54+ }55+ for _, line := range lines {56+ for at := 0; at+1 < len(line); at++ {57+ a, b := line[at], line[at+1]58+ if len(a) < 2 || len(b) < 2 {59+ flush()60+ continue61+ }62+ enter, exit, held := clipSegment(w, a[0], a[1], b[0], b[1])63+ if !held {64+ flush()65+ continue66+ }67+ if len(part) == 0 || part[len(part)-1][0] != enter[0] || part[len(part)-1][1] != enter[1] {68+ flush()69+ part = [][]float64{enter}70+ }71+ part = append(part, exit)72+ if exit[0] != Round7(b[0]) || exit[1] != Round7(b[1]) {73+ flush()74+ }75+ }76+ flush()77+ }78+ return out79+}80+81+// edge is one of the window's four half-planes.82+type edge struct {83+ axis int // 0 clips longitude, 1 latitude84+ bound float6485+ keep float64 // +1 keeps values >= bound, -1 keeps values <= bound86+}87+88+func (w Window) edges() [4]edge {89+ return [4]edge{90+ {0, w.West, 1},91+ {0, w.East, -1},92+ {1, w.South, 1},93+ {1, w.North, -1},94+ }95+}96+97+func (e edge) inside(p []float64) bool {98+ return (p[e.axis]-e.bound)*e.keep >= 099+}100+101+// cross lands on the edge along the segment from a to b, rounded the way102+// every captured coordinate is.103+func (e edge) cross(a, b []float64) []float64 {104+ t := (e.bound - a[e.axis]) / (b[e.axis] - a[e.axis])105+ p := []float64{106+ Round7(a[0] + t*(b[0]-a[0])),107+ Round7(a[1] + t*(b[1]-a[1])),108+ }109+ p[e.axis] = e.bound110+ return p111+}112+113+// clipRing is Sutherland–Hodgman against the window's four edges. The ring114+// arrives closed and leaves closed; nil means nothing of it remains.115+func clipRing(w Window, ring [][]float64) [][]float64 {116+ open := ring117+ if len(open) > 1 && open[0][0] == open[len(open)-1][0] && open[0][1] == open[len(open)-1][1] {118+ open = open[:len(open)-1]119+ }120+ current := make([][]float64, 0, len(open))121+ for _, p := range open {122+ if len(p) >= 2 {123+ current = append(current, p)124+ }125+ }126+ for _, e := range w.edges() {127+ if len(current) == 0 {128+ return nil129+ }130+ next := make([][]float64, 0, len(current)+4)131+ for at, point := range current {132+ previous := current[(at+len(current)-1)%len(current)]133+ switch {134+ case e.inside(point) && e.inside(previous):135+ next = append(next, point)136+ case e.inside(point):137+ next = append(next, e.cross(previous, point), point)138+ case e.inside(previous):139+ next = append(next, e.cross(previous, point))140+ }141+ }142+ current = dedupe(next)143+ }144+ if len(current) < 3 {145+ return nil146+ }147+ return append(current, current[0])148+}149+150+// dedupe drops consecutive repeats, which clipping at a corner introduces151+// and rounding can too.152+func dedupe(points [][]float64) [][]float64 {153+ out := points[:0]154+ for _, p := range points {155+ if len(out) > 0 && out[len(out)-1][0] == p[0] && out[len(out)-1][1] == p[1] {156+ continue157+ }158+ out = append(out, p)159+ }160+ if len(out) > 1 && out[0][0] == out[len(out)-1][0] && out[0][1] == out[len(out)-1][1] {161+ out = out[:len(out)-1]162+ }163+ return out164+}165+166+// clipSegment cuts one segment to the window, answering the kept piece's two167+// ends rounded, or nothing when the segment misses the window entirely.168+func clipSegment(w Window, ax, ay, bx, by float64) (enter, exit []float64, held bool) {169+ t0, t1 := 0.0, 1.0170+ clip := func(delta, room float64) bool {171+ if delta == 0 {172+ return room >= 0173+ }174+ t := room / delta175+ if delta < 0 {176+ if t > t1 {177+ return false178+ }179+ if t > t0 {180+ t0 = t181+ }182+ return true183+ }184+ if t < t0 {185+ return false186+ }187+ if t < t1 {188+ t1 = t189+ }190+ return true191+ }192+ dx, dy := bx-ax, by-ay193+ if !clip(-dx, ax-w.West) || !clip(dx, w.East-ax) ||194+ !clip(-dy, ay-w.South) || !clip(dy, w.North-ay) {195+ return nil, nil, false196+ }197+ enter = []float64{Round7(ax + t0*dx), Round7(ay + t0*dy)}198+ exit = []float64{Round7(ax + t1*dx), Round7(ay + t1*dy)}199+ if enter[0] == exit[0] && enter[1] == exit[1] {200+ return nil, nil, false201+ }202+ return enter, exit, true203+}@@ -0,0 +1,144 @@
1+package arcgismap2+3+import (4+ "encoding/json"5+ "reflect"6+ "testing"7+)8+9+// clipWindow is a plain unit-square window: edges on whole degrees, so every10+// expected crossing is exact.11+var clipWindow = Window{West: 0, South: 0, East: 1, North: 1}12+13+func ring(positions ...[2]float64) [][]float64 {14+ out := make([][]float64, 0, len(positions)+1)15+ for _, p := range positions {16+ out = append(out, []float64{p[0], p[1]})17+ }18+ return append(out, []float64{positions[0][0], positions[0][1]})19+}20+21+func TestClipRingsKeepsAnInsideRing(t *testing.T) {22+ inside := [][][][]float64{{ring(23+ [2]float64{0.2, 0.2}, [2]float64{0.8, 0.2}, [2]float64{0.8, 0.8}, [2]float64{0.2, 0.8},24+ )}}25+ got := ClipRings(clipWindow, inside)26+ if !reflect.DeepEqual(got, inside) {27+ t.Fatalf("an inside ring changed: %v", got)28+ }29+}30+31+func TestClipRingsCutsACrossingRing(t *testing.T) {32+ crossing := [][][][]float64{{ring(33+ [2]float64{0.5, 0.5}, [2]float64{1.5, 0.5}, [2]float64{1.5, 0.7}, [2]float64{0.5, 0.7},34+ )}}35+ got := ClipRings(clipWindow, crossing)36+ if len(got) != 1 || len(got[0]) != 1 {37+ t.Fatalf("crossing ring clipped to %v", got)38+ }39+ want := ring([2]float64{0.5, 0.5}, [2]float64{1, 0.5}, [2]float64{1, 0.7}, [2]float64{0.5, 0.7})40+ if !reflect.DeepEqual(got[0][0], want) {41+ t.Fatalf("crossing ring clipped to %v, want %v", got[0][0], want)42+ }43+}44+45+func TestClipRingsDropsAnOutsideRing(t *testing.T) {46+ outside := [][][][]float64{{ring(47+ [2]float64{2, 2}, [2]float64{3, 2}, [2]float64{3, 3}, [2]float64{2, 3},48+ )}}49+ if got := ClipRings(clipWindow, outside); got != nil {50+ t.Fatalf("an outside ring survived: %v", got)51+ }52+}53+54+func TestClipRingsKeepsAHole(t *testing.T) {55+ holed := [][][][]float64{{56+ ring([2]float64{0.1, 0.1}, [2]float64{0.9, 0.1}, [2]float64{0.9, 0.9}, [2]float64{0.1, 0.9}),57+ ring([2]float64{0.4, 0.4}, [2]float64{0.6, 0.4}, [2]float64{0.6, 0.6}, [2]float64{0.4, 0.6}),58+ }}59+ got := ClipRings(clipWindow, holed)60+ if len(got) != 1 || len(got[0]) != 2 {61+ t.Fatalf("the hole did not survive: %v", got)62+ }63+}64+65+func TestClipRingsDropsAPolygonWhoseGroundIsGone(t *testing.T) {66+ // The outer ring misses the window even though its hole would clip.67+ gone := [][][][]float64{{68+ ring([2]float64{2, 2}, [2]float64{3, 2}, [2]float64{3, 3}, [2]float64{2, 3}),69+ ring([2]float64{0.4, 0.4}, [2]float64{0.6, 0.4}, [2]float64{0.6, 0.6}, [2]float64{0.4, 0.6}),70+ }}71+ if got := ClipRings(clipWindow, gone); got != nil {72+ t.Fatalf("a polygon without ground survived: %v", got)73+ }74+}75+76+func TestClipRingsCoversTheWindowCorner(t *testing.T) {77+ // A huge ring swallowing the whole window clips to the window itself.78+ swallowing := [][][][]float64{{ring(79+ [2]float64{-5, -5}, [2]float64{5, -5}, [2]float64{5, 5}, [2]float64{-5, 5},80+ )}}81+ got := ClipRings(clipWindow, swallowing)82+ if len(got) != 1 || len(got[0]) != 1 || len(got[0][0]) != 5 {83+ t.Fatalf("the window itself should remain, got %v", got)84+ }85+}86+87+func TestClipLinesSplitsAtTheWindow(t *testing.T) {88+ // In at the west, out at the east, back in through the north: two parts.89+ zigzag := [][][]float64{{90+ {-0.5, 0.5}, {0.5, 0.5}, {1.5, 0.5}, {1.5, 2}, {0.5, 2}, {0.5, 0.9},91+ }}92+ got := ClipLines(clipWindow, zigzag)93+ want := [][][]float64{94+ {{0, 0.5}, {0.5, 0.5}, {1, 0.5}},95+ {{0.5, 1}, {0.5, 0.9}},96+ }97+ if !reflect.DeepEqual(got, want) {98+ t.Fatalf("zigzag clipped to %v, want %v", got, want)99+ }100+}101+102+func TestClipLinesKeepsAnInsideLine(t *testing.T) {103+ inside := [][][]float64{{{0.1, 0.1}, {0.5, 0.5}, {0.9, 0.2}}}104+ got := ClipLines(clipWindow, inside)105+ if !reflect.DeepEqual(got, inside) {106+ t.Fatalf("an inside line changed: %v", got)107+ }108+}109+110+func TestClipLinesDropsAnOutsideLine(t *testing.T) {111+ outside := [][][]float64{{{2, 2}, {3, 3}}}112+ if got := ClipLines(clipWindow, outside); got != nil {113+ t.Fatalf("an outside line survived: %v", got)114+ }115+}116+117+// Clipping the same shape twice spells the same bytes, which is what capture118+// determinism asks of every transformation.119+func TestClipIsDeterministic(t *testing.T) {120+ shape := [][][][]float64{{ring(121+ [2]float64{0.123456789, -0.4}, [2]float64{1.7, 0.3}, [2]float64{0.6, 1.9},122+ )}}123+ first, err := json.Marshal(ClipRings(clipWindow, shape))124+ if err != nil {125+ t.Fatal(err)126+ }127+ second, err := json.Marshal(ClipRings(clipWindow, [][][][]float64{{ring(128+ [2]float64{0.123456789, -0.4}, [2]float64{1.7, 0.3}, [2]float64{0.6, 1.9},129+ )}}))130+ if err != nil {131+ t.Fatal(err)132+ }133+ if string(first) != string(second) {134+ t.Fatalf("clipping is not deterministic: %s vs %s", first, second)135+ }136+ again := ClipRings(clipWindow, ClipRings(clipWindow, shape))137+ raw, err := json.Marshal(again)138+ if err != nil {139+ t.Fatal(err)140+ }141+ if string(raw) != string(first) {142+ t.Fatalf("clipping is not idempotent: %s vs %s", raw, first)143+ }144+}@@ -0,0 +1,240 @@
1+package arcgismap2+3+import (4+ "fmt"5+ "math"6+7+ "github.com/FelineStateMachine/atlas/internal/mgdoc"8+ "github.com/FelineStateMachine/atlas/internal/semconv"9+)10+11+// The membership join. A capture that carries the national subwatersheds12+// can say, of each zone the city itself curated, which subwatershed its13+// ground lies in -- computed at translate time from the captured polygons,14+// so re-curating the sentence never needs a re-fetch and the capture stays15+// the one unit of versioning. The claim is made only when all of the16+// zone's sampled ground agrees on one answer: a zone that straddles17+// subwatersheds says nothing rather than something misleading, which is18+// the lesson the excised zoning-rules join left behind -- one earned19+// sentence beats a wall of uncorrelated prose.20+21+// hydroUnit is one subwatershed as captured: its code, its name, and its22+// ground in true degrees.23+type hydroUnit struct {24+ code string25+ name string26+ rings [][][][]float6427+}28+29+// hydroIndex answers which subwatershed a point lies in.30+type hydroIndex struct {31+ units []hydroUnit32+}33+34+// buildHydroIndex reads the captured subwatersheds, if the capture carries35+// them; nil means the capture predates the national enrichment and no zone36+// makes any claim.37+func buildHydroIndex(capture *Capture) *hydroIndex {38+ for _, dataset := range capture.Datasets {39+ if dataset.Slug != SlugSubwatersheds {40+ continue41+ }42+ index := &hydroIndex{}43+ for _, feature := range dataset.Features {44+ code, name := feature.Fields["huc12"], feature.Fields["name"]45+ if code == "" || name == "" || feature.Geometry.Type != GeometryRings {46+ continue47+ }48+ index.units = append(index.units, hydroUnit{49+ code: code, name: name, rings: feature.Geometry.Rings,50+ })51+ }52+ return index53+ }54+ return nil55+}56+57+// locate answers the subwatershed holding a point, if any does. Units58+// arrive normalized, so ties cannot happen except on shared boundaries,59+// where the first unit in capture order answers deterministically.60+func (h *hydroIndex) locate(lon, lat float64) (hydroUnit, bool) {61+ if h == nil {62+ return hydroUnit{}, false63+ }64+ for _, unit := range h.units {65+ if pointInRings(lon, lat, unit.rings) {66+ return unit, true67+ }68+ }69+ return hydroUnit{}, false70+}71+72+// pointInRings is even-odd ray casting over MultiPolygon nesting: crossing73+// an outer ring enters ground, crossing a hole leaves it again.74+func pointInRings(lon, lat float64, polygons [][][][]float64) bool {75+ inside := false76+ for _, polygon := range polygons {77+ for _, ring := range polygon {78+ for at := 0; at+1 < len(ring); at++ {79+ a, b := ring[at], ring[at+1]80+ if len(a) < 2 || len(b) < 2 {81+ continue82+ }83+ if (a[1] > lat) == (b[1] > lat) {84+ continue85+ }86+ crossing := a[0] + (lat-a[1])/(b[1]-a[1])*(b[0]-a[0])87+ if lon < crossing {88+ inside = !inside89+ }90+ }91+ }92+ }93+ return inside94+}95+96+// featureAnchor is one feature's representative point in true degrees: the97+// centroid of its largest outer ring, or the middle vertex of its longest98+// line. A point feature is its own anchor, though no point dataset zones.99+func featureAnchor(g Geometry) (lon, lat float64, ok bool) {100+ switch g.Type {101+ case GeometryPoint:102+ if len(g.Point) < 2 {103+ return 0, 0, false104+ }105+ return g.Point[0], g.Point[1], true106+ case GeometryRings:107+ best := 0.0108+ for _, polygon := range g.Rings {109+ if len(polygon) == 0 {110+ continue111+ }112+ x, y, area, held := ringCentroid(polygon[0])113+ if held && area > best {114+ best, lon, lat, ok = area, x, y, true115+ }116+ }117+ return lon, lat, ok118+ case GeometryLines:119+ var longest []float64120+ best := 0.0121+ for _, line := range g.Lines {122+ length := 0.0123+ for at := 0; at+1 < len(line); at++ {124+ if len(line[at]) < 2 || len(line[at+1]) < 2 {125+ continue126+ }127+ length += math.Hypot(line[at+1][0]-line[at][0], line[at+1][1]-line[at][1])128+ }129+ if length > best && len(line) > 0 {130+ best, longest = length, line[len(line)/2]131+ }132+ }133+ if len(longest) < 2 {134+ return 0, 0, false135+ }136+ return longest[0], longest[1], true137+ }138+ return 0, 0, false139+}140+141+// ringCentroid is the shoelace centroid of one ring, with the enclosed142+// area to rank rings by.143+func ringCentroid(ring [][]float64) (lon, lat, area float64, ok bool) {144+ var doubled, sumX, sumY float64145+ for at := 0; at+1 < len(ring); at++ {146+ a, b := ring[at], ring[at+1]147+ if len(a) < 2 || len(b) < 2 {148+ return 0, 0, 0, false149+ }150+ cross := a[0]*b[1] - b[0]*a[1]151+ doubled += cross152+ sumX += (a[0] + b[0]) * cross153+ sumY += (a[1] + b[1]) * cross154+ }155+ if doubled == 0 {156+ return 0, 0, 0, false157+ }158+ return sumX / (3 * doubled), sumY / (3 * doubled), math.Abs(doubled) / 2, true159+}160+161+// hydroClaims tracks, per zone, the subwatersheds its features' ground162+// lands in, and whether any of it landed nowhere.163+type hydroClaims struct {164+ codes map[string]bool165+ missed bool166+}167+168+func (c *hydroClaims) observe(index *hydroIndex, g Geometry) {169+ for _, position := range claimPositions(g) {170+ unit, found := index.locate(position[0], position[1])171+ if !found {172+ c.missed = true173+ continue174+ }175+ if c.codes == nil {176+ c.codes = map[string]bool{}177+ }178+ c.codes[unit.code] = true179+ }180+}181+182+// claimSamples is how many boundary positions a feature volunteers beside183+// its anchor. The anchor alone would let a zone made of one citywide184+// polygon claim whichever subwatershed its centroid happens to sit in;185+// sampling the boundary catches the spread, and the cost stays flat no186+// matter how detailed the ring.187+const claimSamples = 16188+189+// claimPositions is the ground a feature answers for: its anchor, and an190+// even sample of its outer-ring or line vertices.191+func claimPositions(g Geometry) [][2]float64 {192+ var out [][2]float64193+ if lon, lat, ok := featureAnchor(g); ok {194+ out = append(out, [2]float64{lon, lat})195+ }196+ var boundary [][]float64197+ switch g.Type {198+ case GeometryRings:199+ for _, polygon := range g.Rings {200+ if len(polygon) > 0 {201+ boundary = append(boundary, polygon[0]...)202+ }203+ }204+ case GeometryLines:205+ for _, line := range g.Lines {206+ boundary = append(boundary, line...)207+ }208+ }209+ step := max((len(boundary)+claimSamples-1)/claimSamples, 1)210+ for at := 0; at < len(boundary); at += step {211+ if len(boundary[at]) >= 2 {212+ out = append(out, [2]float64{boundary[at][0], boundary[at][1]})213+ }214+ }215+ return out216+}217+218+// apply writes the membership onto a zone when every feature agreed on one219+// subwatershed: a short sentence for the card, and the code as the220+// machine-readable key beside it.221+func (c *hydroClaims) apply(index *hydroIndex, zone *mgdoc.Region) {222+ if c.missed || len(c.codes) != 1 {223+ return224+ }225+ var code string226+ for held := range c.codes {227+ code = held228+ }229+ for _, unit := range index.units {230+ if unit.code != code {231+ continue232+ }233+ zone.Description = scrub(fmt.Sprintf("Lies in the %s subwatershed (HUC %s).", unit.name, unit.code))234+ if zone.Attrs == nil {235+ zone.Attrs = map[string]string{}236+ }237+ zone.Attrs[semconv.KeyHydroHUC12] = unit.code238+ return239+ }240+}@@ -0,0 +1,106 @@
1+package arcgismap2+3+import (4+ "math"5+ "testing"6+)7+8+func squareRings(west, south, size float64) [][][][]float64 {9+ return [][][][]float64{{{10+ {west, south}, {west + size, south}, {west + size, south + size}, {west, south + size}, {west, south},11+ }}}12+}13+14+func TestPointInRings(t *testing.T) {15+ holed := [][][][]float64{{16+ {{0, 0}, {10, 0}, {10, 10}, {0, 10}, {0, 0}},17+ {{4, 4}, {6, 4}, {6, 6}, {4, 6}, {4, 4}},18+ }}19+ cases := []struct {20+ lon, lat float6421+ want bool22+ }{23+ {5, 2, true}, // on the ground24+ {5, 5, false}, // in the hole25+ {15, 5, false}, // outside26+ {-1, -1, false},27+ }28+ for _, tc := range cases {29+ if got := pointInRings(tc.lon, tc.lat, holed); got != tc.want {30+ t.Fatalf("pointInRings(%v,%v) = %v", tc.lon, tc.lat, got)31+ }32+ }33+}34+35+func TestFeatureAnchor(t *testing.T) {36+ lon, lat, ok := featureAnchor(Geometry{Type: GeometryRings, Rings: squareRings(2, 4, 2)})37+ if !ok || math.Abs(lon-3) > 1e-9 || math.Abs(lat-5) > 1e-9 {38+ t.Fatalf("square anchors at (%v,%v), %v", lon, lat, ok)39+ }40+ // Two polygons: the larger one's centroid answers.41+ both := append(squareRings(0, 0, 1), squareRings(10, 10, 3)...)42+ lon, lat, ok = featureAnchor(Geometry{Type: GeometryRings, Rings: both})43+ if !ok || math.Abs(lon-11.5) > 1e-9 || math.Abs(lat-11.5) > 1e-9 {44+ t.Fatalf("the larger polygon should anchor, got (%v,%v)", lon, lat)45+ }46+ // The longest line's middle vertex answers.47+ lon, lat, ok = featureAnchor(Geometry{Type: GeometryLines, Lines: [][][]float64{48+ {{0, 0}, {0.1, 0}},49+ {{5, 5}, {6, 5}, {7, 5}},50+ }})51+ if !ok || lon != 6 || lat != 5 {52+ t.Fatalf("line anchors at (%v,%v), %v", lon, lat, ok)53+ }54+ if _, _, ok := featureAnchor(Geometry{Type: GeometryRings}); ok {55+ t.Fatal("an empty geometry anchored")56+ }57+}58+59+func TestLocateFindsTheHoldingUnit(t *testing.T) {60+ index := &hydroIndex{units: []hydroUnit{61+ {code: "101900030301", name: "Walnut Creek", rings: squareRings(0, 0, 1)},62+ {code: "101900030304", name: "Big Dry Creek", rings: squareRings(1, 0, 1)},63+ }}64+ unit, found := index.locate(1.5, 0.5)65+ if !found || unit.name != "Big Dry Creek" {66+ t.Fatalf("located %+v, %v", unit, found)67+ }68+ if _, found := index.locate(5, 5); found {69+ t.Fatal("a point in no unit located one")70+ }71+ var none *hydroIndex72+ if _, found := none.locate(0, 0); found {73+ t.Fatal("a nil index located a unit")74+ }75+}76+77+func TestClaimPositionsSamplesTheBoundary(t *testing.T) {78+ // A long line: the anchor plus an even sample, never the whole ring.79+ line := make([][]float64, 100)80+ for at := range line {81+ line[at] = []float64{float64(at), 0}82+ }83+ positions := claimPositions(Geometry{Type: GeometryLines, Lines: [][][]float64{line}})84+ if len(positions) < 2 || len(positions) > claimSamples+2 {85+ t.Fatalf("%d positions sampled", len(positions))86+ }87+}88+89+func TestBuildHydroIndexReadsTheCapturedSubwatersheds(t *testing.T) {90+ capture := fixtureCapture()91+ index := buildHydroIndex(&capture)92+ if index == nil || len(index.units) != 2 {93+ t.Fatalf("index is %+v", index)94+ }95+ without := fixtureCapture()96+ kept := without.Datasets[:0]97+ for _, dataset := range without.Datasets {98+ if dataset.Slug != SlugSubwatersheds {99+ kept = append(kept, dataset)100+ }101+ }102+ without.Datasets = kept103+ if buildHydroIndex(&without) != nil {104+ t.Fatal("a capture without subwatersheds built an index")105+ }106+}@@ -0,0 +1,108 @@
1+package arcgismap2+3+import (4+ "strconv"5+ "strings"6+)7+8+// The national enrichment. A state GIS viewer is usually a branded window9+// onto federal datasets -- Colorado's water viewer credits the USGS10+// Watershed Boundary Dataset and National Hydrography Dataset by name -- so11+// the generic move is to ask the federal services directly, by any city's12+// own bounding box, and every curated city gains the same layers without a13+// line of per-city curation.14+//15+// Both services are plain ArcGIS MapServer query endpoints on16+// hydro.nationalmap.gov. Field spellings differ per layer -- the WBD layers17+// speak lowercase, the NHD waterbody layer uppercase -- and were verified18+// against the live services, the same discipline every city table keeps.19+20+const (21+ wbdServer = "https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer"22+ nhdServer = "https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer"23+24+ // SlugSubwatersheds names the captured dataset the membership join25+ // reads, so the join and the table cannot drift apart on a spelling.26+ SlugSubwatersheds = "subwatersheds"27+)28+29+// National is the enrichment every curated city receives: the watersheds30+// its ground drains through, in two grains of the USGS hydrologic unit31+// hierarchy, and the named streams and waterbodies that do the draining.32+// A padded city window intersects a handful of hydrologic units and tens33+// of named streams -- national extent, city-sized answers.34+var National = []Dataset{35+ {36+ Slug: "watersheds", Title: "Watersheds",37+ Server: wbdServer, Layer: 5,38+ Geometry: "polygon",39+ Keep: []string{"objectid", "huc10", "name"},40+ IDOf: hucID("huc10"),41+ ZoneOf: hucZone("huc10", "Watershed"),42+ },43+ {44+ Slug: SlugSubwatersheds, Title: "Subwatersheds",45+ Server: wbdServer, Layer: 6,46+ Geometry: "polygon",47+ Keep: []string{"objectid", "huc12", "name"},48+ IDOf: hucID("huc12"),49+ ZoneOf: hucZone("huc12", "Subwatershed"),50+ },51+ {52+ Slug: "streams", Title: "Streams",53+ Server: nhdServer, Layer: 6,54+ Geometry: "line",55+ Keep: []string{"OBJECTID", "gnis_name"},56+ // Named flowlines only: the layer also carries every unnamed ditch57+ // and culvert, which is noise to a reader and no loss to the map.58+ Where: "gnis_name IS NOT NULL",59+ ZoneOf: func(f Fields) ZoneKey {60+ name := strings.TrimSpace(f["gnis_name"])61+ if name == "" {62+ return ZoneKey{}63+ }64+ return ZoneKey{Key: slugify(name), Title: name, Subtitle: "Stream"}65+ },66+ StrokeWidth: 10,67+ Role: "water",68+ },69+ {70+ Slug: "waterbodies", Title: "Waterbodies",71+ Server: nhdServer, Layer: 12,72+ Geometry: "polygon",73+ Keep: []string{"OBJECTID", "GNIS_NAME"},74+ // Every pond draws into the basemap; only the named ones earn a75+ // zone anyone could search for.76+ ZoneOf: func(f Fields) ZoneKey {77+ name := strings.TrimSpace(f["GNIS_NAME"])78+ if name == "" {79+ return ZoneKey{}80+ }81+ return ZoneKey{Key: slugify(name), Title: name, Subtitle: "Waterbody"}82+ },83+ Role: "water",84+ },85+}86+87+// hucID reads the hydrologic unit code as the feature's identity: the code88+// is nationally unique and survives upstream refreshes, where the layer's89+// object ids are load artifacts that churn. A leading zero -- regions one90+// through nine -- drops in the number without costing uniqueness.91+func hucID(field string) func(Fields) (int64, bool) {92+ return func(f Fields) (int64, bool) {93+ code, err := strconv.ParseInt(f[field], 10, 64)94+ return code, err == nil && code > 095+ }96+}97+98+// hucZone folds a WBD layer into one zone per hydrologic unit, titled by99+// the unit's own name with the code kept legible beside the grain.100+func hucZone(field, grain string) func(Fields) ZoneKey {101+ return func(f Fields) ZoneKey {102+ code, name := strings.TrimSpace(f[field]), strings.TrimSpace(f["name"])103+ if code == "" || name == "" {104+ return ZoneKey{}105+ }106+ return ZoneKey{Key: code, Title: name, Subtitle: grain + " · HUC " + code}107+ }108+}@@ -29,6 +29,7 @@ breaks; new keys arrive as `experimental` and earn `stable`.
2929 | `atlas.geometry.radius_km` | world | decimal string | experimental | The body's mean radius. |3030 | `atlas.geo.lat` | location | decimal degrees | experimental | True planetary latitude as published (planetocentric). Provenance and card material; rendering derives from the map-level mapping. |3131 | `atlas.geo.lon` | location | decimal degrees | experimental | True planetary longitude as published (east-positive). |32+| `atlas.hydro.huc12` | zone | twelve digits | experimental | The USGS twelve-digit hydrologic unit — the subwatershed — the zone's ground lies wholly within, from the national Watershed Boundary Dataset. A zone spanning subwatersheds carries no key rather than a misleading one. |3233 | `atlas.stroke.width_px` | zone | positive decimal (world px) | experimental | The ground width of a zone whose features are lines: a trail is a line and a weight, and the weight lets a reader draw the path as one continuous stroke. |3334 | `atlas.category.key` | category | slug | experimental | A category's merge identity across sources. Absent, the icon key stands in. |3435 @@ -116,6 +116,13 @@ const (
116116 KeyGeoLat = "atlas.geo.lat"117117 KeyGeoLon = "atlas.geo.lon"118118 119+ // HydroHUC12 names the USGS twelve-digit hydrologic unit -- the120+ // subwatershed -- a zone's ground lies wholly within, from the national121+ // Watershed Boundary Dataset, so watershed membership is machine-readable122+ // beside whatever sentence the zone's card shows. A zone spanning123+ // subwatersheds carries no key rather than a misleading one.124+ KeyHydroHUC12 = "atlas.hydro.huc12"125+119126 // StrokeWidthPx is the ground width of a zone whose features are lines120127 // rather than areas, in world pixels: a trail is a line and a weight,121128 // and declaring the weight lets a reader draw the path as one@@ -196,6 +203,18 @@ func positiveDecimal(value string) error {
196203 return nil197204 }198205 206+func huc12(value string) error {207+ if len(value) != 12 {208+ return fmt.Errorf("%q is not twelve digits", value)209+ }210+ for _, r := range value {211+ if r < '0' || r > '9' {212+ return fmt.Errorf("%q is not twelve digits", value)213+ }214+ }215+ return nil216+}217+199218 func setName(value string) error {200219 set, name, found := strings.Cut(value, "/")201220 if !found || set == "" || name == "" {@@ -239,6 +258,7 @@ var registry = map[string]definition{
239258 KeyGeometryRadiusKM: {EntityWorld, Experimental, decimal},240259 KeyGeoLat: {EntityLocation, Experimental, decimal},241260 KeyGeoLon: {EntityLocation, Experimental, decimal},261+ KeyHydroHUC12: {EntityZone, Experimental, huc12},242262 KeyStrokeWidthPx: {EntityZone, Experimental, positiveDecimal},243263 KeyCategoryKey: {EntityCategory, Experimental, slug},244264 }@@ -131,7 +131,7 @@ func TestRegistryAgreesWithItsDocument(t *testing.T) {
131131 if err != nil {132132 t.Fatal(err)133133 }134- row := regexp.MustCompile("(?m)^\\| `(atlas\\.[a-z_.]+)` \\| (bundle|world|zone|category|location) \\|.*\\| (stable|experimental) \\|")134+ row := regexp.MustCompile("(?m)^\\| `(atlas\\.[a-z0-9_.]+)` \\| (bundle|world|zone|category|location) \\|.*\\| (stable|experimental) \\|")135135 documented := make(map[string][2]string)136136 for _, match := range row.FindAllStringSubmatch(string(doc), -1) {137137 documented[match[1]] = [2]string{match[2], match[3]}@@ -7,7 +7,11 @@
77 // This is the curated-cities importer, and says so. Each city in the88 // arcgismap table names verified dataset identities, field names, and a99 // bounding box; an unknown city is refused at the door rather than captured10-// through unverified curation.10+// through unverified curation. Beside the city's hub, every crawl asks the11+// USGS national services for the city's window -- the national enrichment12+// in arcgismap's National table, fetched through the MapServer query path13+// in mapserver.go -- so the capture holds the city and its watershed14+// context as one document.1115 //1216 // Versioning is the mode's whole reason. Each crawl day registers its own1317 // map directory, so the archive keeps every day the data moved and the@@ -73,10 +77,17 @@ func runArcgis(ctx context.Context, fetcher *fetcher, o options) error {
7377 Basemap: arcgismap.MapConfig{MaxZoom: maxZoom, Extension: "png"},7478 }7579 76- fmt.Printf("\n== ArcGIS Hub %s / %s (%d datasets)\n", city.Title, date, len(city.Datasets))80+ datasets := city.AllDatasets()81+ fmt.Printf("\n== ArcGIS Hub %s / %s (%d datasets)\n", city.Title, date, len(datasets))7782 pinCount := 078- for _, dataset := range city.Datasets {79- features, err := fetchArcgisDataset(ctx, fetcher, city, dataset)83+ for _, dataset := range datasets {84+ var features []arcgismap.Feature85+ var err error86+ if dataset.Server != "" {87+ features, err = fetchMapServerDataset(ctx, fetcher, window, dataset)88+ } else {89+ features, err = fetchArcgisDataset(ctx, fetcher, city, dataset)90+ }8091 if err != nil {8192 return fmt.Errorf("%s: %w", dataset.Slug, err)8293 }@@ -514,9 +525,10 @@ func renderArcgisBasemap(
514525 // capture's own window, so the raster and the pins can never disagree about515526 // where the ground is.516527 func basemapFeatures(city arcgismap.City, capture *arcgismap.Capture) []basemap.Feature {517- curated := make(map[string]*arcgismap.Dataset, len(city.Datasets))518- for at := range city.Datasets {519- curated[city.Datasets[at].Slug] = &city.Datasets[at]528+ datasets := city.AllDatasets()529+ curated := make(map[string]*arcgismap.Dataset, len(datasets))530+ for at := range datasets {531+ curated[datasets[at].Slug] = &datasets[at]520532 }521533 var features []basemap.Feature522534 for _, dataset := range capture.Datasets {@@ -0,0 +1,193 @@
1+// The national fetch path. Hub datasets arrive through a download API that2+// stages a whole file server-side; a national layer on a plain ArcGIS3+// MapServer offers no such staging, so it is paged through the query4+// endpoint instead -- the city's padded window as the envelope, true degrees5+// asked for by name, pages ordered by the layer's own row identity until6+// one comes back short. Fetched geometry is cut to the window at capture7+// time: a watershed intersecting the window may extend a hundred windows8+// past it, and the trimmed shape is the honest city-sized answer.9+package main10+11+import (12+ "context"13+ "encoding/json"14+ "fmt"15+ "net/url"16+ "strconv"17+ "strings"18+19+ "github.com/FelineStateMachine/atlas/internal/arcgismap"20+)21+22+// mapServerPage is the most rows one query asks for, matching the services'23+// own maxRecordCount.24+const mapServerPage = 200025+26+// fetchMapServerDataset takes one national layer through the query endpoint,27+// page by page, and answers features already cut to the window.28+func fetchMapServerDataset(29+ ctx context.Context,30+ fetcher *fetcher,31+ window arcgismap.Window,32+ dataset arcgismap.Dataset,33+) ([]arcgismap.Feature, error) {34+ var features []arcgismap.Feature35+ claimed := map[int64]bool{}36+ for offset := 0; ; offset += mapServerPage {37+ address := mapServerQueryURL(dataset, window, offset)38+ body, _, err := fetcher.get(ctx, address)39+ if err != nil {40+ return nil, err41+ }42+ page, more, err := decodeMapServerPage(dataset, window, body, address)43+ if err != nil {44+ return nil, err45+ }46+ for _, feature := range page {47+ // A row captured twice means the pages shifted under the48+ // crawl; a capture must not hold it two ways, and Normalize's49+ // sort would bury the evidence.50+ if claimed[feature.ID] {51+ return nil, fmt.Errorf("feature %d arrived twice; pages shifted mid-crawl: %s",52+ feature.ID, address)53+ }54+ claimed[feature.ID] = true55+ }56+ features = append(features, page...)57+ if !more {58+ return features, nil59+ }60+ }61+}62+63+// mapServerQueryURL spells one page's request: the window as the envelope,64+// degrees in and degrees out, fields limited to the curation, rows ordered65+// by the layer's row identity so the pages tile the layer exactly once.66+func mapServerQueryURL(dataset arcgismap.Dataset, window arcgismap.Window, offset int) string {67+ where := dataset.Where68+ if where == "" {69+ where = "1=1"70+ }71+ envelope := strings.Join([]string{72+ strconv.FormatFloat(window.West, 'f', -1, 64),73+ strconv.FormatFloat(window.South, 'f', -1, 64),74+ strconv.FormatFloat(window.East, 'f', -1, 64),75+ strconv.FormatFloat(window.North, 'f', -1, 64),76+ }, ",")77+ values := url.Values{78+ "f": {"geojson"},79+ "where": {where},80+ "geometry": {envelope},81+ "geometryType": {"esriGeometryEnvelope"},82+ "inSR": {"4326"},83+ "spatialRel": {"esriSpatialRelIntersects"},84+ "outSR": {"4326"},85+ "outFields": {strings.Join(dataset.Keep, ",")},86+ "geometryPrecision": {"7"},87+ "orderByFields": {dataset.Keep[0]},88+ "resultOffset": {strconv.Itoa(offset)},89+ "resultRecordCount": {strconv.Itoa(mapServerPage)},90+ }91+ return fmt.Sprintf("%s/%d/query?%s", dataset.Server, dataset.Layer, values.Encode())92+}93+94+// decodeMapServerPage reads one page into capture features, cut to the95+// window, and says whether another page waits. ArcGIS answers errors as96+// HTTP 200 with an error member, and spells the truncation flag at the97+// collection's top level or under its properties depending on the vintage.98+func decodeMapServerPage(99+ dataset arcgismap.Dataset,100+ window arcgismap.Window,101+ body []byte,102+ address string,103+) ([]arcgismap.Feature, bool, error) {104+ var page struct {105+ Error *struct {106+ Code int `json:"code"`107+ Message string `json:"message"`108+ } `json:"error"`109+ Type string `json:"type"`110+ ExceededTransferLimit bool `json:"exceededTransferLimit"`111+ Properties struct {112+ ExceededTransferLimit bool `json:"exceededTransferLimit"`113+ } `json:"properties"`114+ Features []struct {115+ Properties map[string]any `json:"properties"`116+ Geometry struct {117+ Type string `json:"type"`118+ Coordinates json.RawMessage `json:"coordinates"`119+ } `json:"geometry"`120+ } `json:"features"`121+ }122+ if err := json.Unmarshal(body, &page); err != nil {123+ return nil, false, fmt.Errorf("decode %s: %w", address, err)124+ }125+ if page.Error != nil {126+ return nil, false, fmt.Errorf("%s answered error %d: %s",127+ address, page.Error.Code, page.Error.Message)128+ }129+ if page.Type != "FeatureCollection" {130+ return nil, false, fmt.Errorf("%s is a %q, not a FeatureCollection", address, page.Type)131+ }132+133+ features := make([]arcgismap.Feature, 0, len(page.Features))134+ for at, raw := range page.Features {135+ fields := arcgismap.Fields{}136+ for _, keep := range dataset.Keep {137+ if value := arcgismap.FieldString(raw.Properties[keep]); value != "" {138+ fields[keep] = value139+ }140+ }141+ id, err := mapServerFeatureID(dataset, fields, raw.Properties)142+ if err != nil {143+ return nil, false, fmt.Errorf("feature %d of %s: %w", at, address, err)144+ }145+ geometry, err := arcgisGeometry(dataset.Geometry, raw.Geometry.Type, raw.Geometry.Coordinates)146+ if err != nil {147+ return nil, false, fmt.Errorf("feature %d: %w", id, err)148+ }149+ if geometry = clipToWindow(window, geometry); geometry == nil {150+ continue151+ }152+ features = append(features, arcgismap.Feature{ID: id, Fields: fields, Geometry: *geometry})153+ }154+ more := page.ExceededTransferLimit || page.Properties.ExceededTransferLimit ||155+ len(page.Features) == mapServerPage156+ return features, more, nil157+}158+159+// mapServerFeatureID is the capture's identity for one national row: the160+// curated derivation where the table declares one, the layer's own object161+// id otherwise.162+func mapServerFeatureID(dataset arcgismap.Dataset, fields arcgismap.Fields, properties map[string]any) (int64, error) {163+ if dataset.IDOf == nil {164+ return arcgisObjectID(properties)165+ }166+ id, ok := dataset.IDOf(fields)167+ if !ok {168+ return 0, fmt.Errorf("no identity in %v", fields)169+ }170+ return id, nil171+}172+173+// clipToWindow cuts a normalized geometry to the window, answering nil when174+// nothing of it remains. Points pass untouched: no national layer captures175+// them yet, and a stray point is the stray tolerance's business, not this.176+func clipToWindow(window arcgismap.Window, geometry *arcgismap.Geometry) *arcgismap.Geometry {177+ if geometry == nil {178+ return nil179+ }180+ switch geometry.Type {181+ case arcgismap.GeometryRings:182+ geometry.Rings = arcgismap.ClipRings(window, geometry.Rings)183+ if len(geometry.Rings) == 0 {184+ return nil185+ }186+ case arcgismap.GeometryLines:187+ geometry.Lines = arcgismap.ClipLines(window, geometry.Lines)188+ if len(geometry.Lines) == 0 {189+ return nil190+ }191+ }192+ return geometry193+}@@ -0,0 +1,120 @@
1+package main2+3+import (4+ "net/url"5+ "strings"6+ "testing"7+8+ "github.com/FelineStateMachine/atlas/internal/arcgismap"9+)10+11+var testWindow = arcgismap.Window{West: -105.2, South: 39.8, East: -104.9, North: 40.0}12+13+func testDataset() arcgismap.Dataset {14+ for _, dataset := range arcgismap.National {15+ if dataset.Slug == "streams" {16+ return dataset17+ }18+ }19+ panic("streams is not curated")20+}21+22+func TestMapServerQueryURLSpellsThePage(t *testing.T) {23+ address := mapServerQueryURL(testDataset(), testWindow, 4000)24+ parsed, err := url.Parse(address)25+ if err != nil {26+ t.Fatal(err)27+ }28+ if !strings.HasSuffix(parsed.Path, "/6/query") {29+ t.Fatalf("path is %q", parsed.Path)30+ }31+ values := parsed.Query()32+ wants := map[string]string{33+ "f": "geojson",34+ "where": "gnis_name IS NOT NULL",35+ "geometry": "-105.2,39.8,-104.9,40",36+ "geometryType": "esriGeometryEnvelope",37+ "inSR": "4326",38+ "outSR": "4326",39+ "spatialRel": "esriSpatialRelIntersects",40+ "outFields": "OBJECTID,gnis_name",41+ "geometryPrecision": "7",42+ "orderByFields": "OBJECTID",43+ "resultOffset": "4000",44+ "resultRecordCount": "2000",45+ }46+ for key, want := range wants {47+ if got := values.Get(key); got != want {48+ t.Fatalf("%s is %q, want %q", key, got, want)49+ }50+ }51+ // The where-clause must actually be escaped on the wire.52+ if strings.Contains(address, "gnis_name IS") {53+ t.Fatalf("where rides unescaped: %s", address)54+ }55+}56+57+func TestDecodeMapServerPageReadsFeaturesAndClips(t *testing.T) {58+ body := []byte(`{"type":"FeatureCollection","features":[59+ {"properties":{"OBJECTID":7,"gnis_name":"Big Dry Creek"},60+ "geometry":{"type":"LineString","coordinates":[[-105.1,39.9],[-105.0,39.9]]}},61+ {"properties":{"OBJECTID":8,"gnis_name":"Elsewhere Creek"},62+ "geometry":{"type":"LineString","coordinates":[[-99.0,39.9],[-98.9,39.9]]}}63+ ]}`)64+ features, more, err := decodeMapServerPage(testDataset(), testWindow, body, "test")65+ if err != nil {66+ t.Fatal(err)67+ }68+ if more {69+ t.Fatal("a short page claims another waits")70+ }71+ if len(features) != 1 || features[0].ID != 7 {72+ t.Fatalf("features are %+v; the out-of-window line should be cut away", features)73+ }74+ if features[0].Fields["gnis_name"] != "Big Dry Creek" {75+ t.Fatalf("fields are %v", features[0].Fields)76+ }77+}78+79+func TestDecodeMapServerPageReadsTheTruncationFlag(t *testing.T) {80+ spellings := []string{81+ `{"type":"FeatureCollection","exceededTransferLimit":true,"features":[]}`,82+ `{"type":"FeatureCollection","properties":{"exceededTransferLimit":true},"features":[]}`,83+ }84+ for _, body := range spellings {85+ _, more, err := decodeMapServerPage(testDataset(), testWindow, []byte(body), "test")86+ if err != nil {87+ t.Fatal(err)88+ }89+ if !more {90+ t.Fatalf("truncation went unread in %s", body)91+ }92+ }93+}94+95+func TestDecodeMapServerPageRefusesAnErrorAnswer(t *testing.T) {96+ body := []byte(`{"error":{"code":400,"message":"Invalid parameters"}}`)97+ if _, _, err := decodeMapServerPage(testDataset(), testWindow, body, "test"); err == nil {98+ t.Fatal("an error answer decoded as a page")99+ }100+}101+102+func TestMapServerFeatureIDPrefersTheCuratedDerivation(t *testing.T) {103+ var subwatersheds arcgismap.Dataset104+ for _, dataset := range arcgismap.National {105+ if dataset.Slug == arcgismap.SlugSubwatersheds {106+ subwatersheds = dataset107+ }108+ }109+ id, err := mapServerFeatureID(subwatersheds,110+ arcgismap.Fields{"huc12": "101900030304"}, map[string]any{"objectid": float64(9)})111+ if err != nil {112+ t.Fatal(err)113+ }114+ if id != 101900030304 {115+ t.Fatalf("id is %d, want the hydrologic unit code", id)116+ }117+ if _, err := mapServerFeatureID(subwatersheds, arcgismap.Fields{}, nil); err == nil {118+ t.Fatal("a row without its code passed")119+ }120+}