The seam the zoning-rules excise left open finally has data that earns it. A capture carrying the national subwatersheds can say, of each zone the city itself curated, which one the zone's ground lies in -- joined at translate time from the captured polygons, so re-curating the sentence never needs a re-fetch and the capture stays the one unit of versioning. The claim is guarded the way that excise taught: a zone speaks only when its anchor and an even sample of its boundary all land in the same subwatershed. A zone that straddles two says nothing; ground that lands in no captured unit says nothing; the national zones never describe themselves. What survives the guard is one earned sentence on the card and the twelve-digit code in atlas.hydro.huc12 beside it.
- commit
- 0ef7f7c8ba848a51920817b743f3af2ad53e69cb
- parent
- 80c32362534153e4e5092a3163922dc39bf2cfc7
- author
- Dami <dami.etoile@gmail.com>
- kind
- feat risk score 1.4 of 3, Jev's estimate
| M | internal/ | 25 | +++- |
| M | internal/ | 73 | +++++++++- |
| A | internal/ | 240 | ++++++++++++++++++++++++++++++ |
| A | internal/ | 106 | +++++++++++++ |
@@ -547,7 +547,7 @@ func Translate(doc []byte) ([]byte, error) {
547547 if err != nil {548548 return nil, err549549 }550- regions, err := buildRegions(&capture, byOrder, ids)550+ regions, err := buildRegions(&capture, byOrder, ids, buildHydroIndex(&capture))551551 if err != nil {552552 return nil, err553553 }@@ -758,8 +758,10 @@ const zoneLimit = 256
758758 // buildRegions folds the polygon and line datasets into zones: every759759 // feature whose ZoneOf names a bucket lands its ground in that bucket's760760 // region, and line features widen into ribbon polygons first, because a761-// zone is ground and a line has none.762-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) {763765 regions := []mgdoc.Region{}764766 for _, pair := range pairs {765767 curated, data := pair.curated, pair.data@@ -767,12 +769,21 @@ func buildRegions(capture *Capture, pairs []pairing, ids *mgdoc.IDSpace) ([]mgdo
767769 continue768770 }769771 buckets := make(map[string]*mgdoc.Region)772+ claims := make(map[string]*hydroClaims)770773 var order []string771774 for _, feature := range data.Features {772775 key := curated.ZoneOf(feature.Fields)773776 if key.Key == "" {774777 continue775778 }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+ }776787 zone, made := buckets[key.Key]777788 if !made {778789 if len(buckets) >= zoneLimit {@@ -807,9 +818,13 @@ func buildRegions(capture *Capture, pairs []pairing, ids *mgdoc.IDSpace) ([]mgdo
807818 }808819 sort.Strings(order)809820 for _, key := range order {810- if len(buckets[key].Features) > 0 {811- 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])812826 }827+ regions = append(regions, *buckets[key])813828 }814829 }815830 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 @@ -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+}