1package arcgishub23import (4 "fmt"5 "math"67 "github.com/FelineStateMachine/atlas/format/semconv"8 "github.com/FelineStateMachine/atlas/internal/generate/doc"9)1011// The membership join. A capture that carries the national subwatersheds can12// say, of each zone the city itself curated, which subwatershed its ground lies13// in -- computed here from the captured polygons, so re-curating the sentence14// never needs a re-fetch and the capture stays the one unit of versioning.15//16// The claim is made only when all of a zone's sampled ground agrees on one17// answer. A zone that straddles subwatersheds says nothing rather than18// something misleading: one earned sentence beats a wall of uncorrelated prose.1920// hydroUnit is one subwatershed as captured: its code, its name, and its ground21// in true degrees.22type hydroUnit struct {23 code string24 name string25 rings [][][][]float6426}2728// hydroIndex answers which subwatershed a position lies in.29type hydroIndex struct {30 units []hydroUnit31}3233// buildHydroIndex reads the captured subwatersheds, if the capture carries34// them; nil means the capture predates the national enrichment and no zone35// makes any claim.36func buildHydroIndex(raw *capture) *hydroIndex {37 for _, captured := range raw.Datasets {38 if captured.Slug != slugSubwatersheds {39 continue40 }41 index := &hydroIndex{}42 for _, row := range captured.Features {43 code, name := row.Fields["huc12"], row.Fields["name"]44 if code == "" || name == "" || row.Geometry.Type != geometryRings {45 continue46 }47 index.units = append(index.units, hydroUnit{code: code, name: name, rings: row.Geometry.Rings})48 }49 return index50 }51 return nil52}5354// locate answers the subwatershed holding a position, if any does. Units arrive55// in capture order, so a position on a shared boundary answers deterministically56// rather than by whichever unit was tried first.57func (h *hydroIndex) locate(lon, lat float64) (hydroUnit, bool) {58 if h == nil {59 return hydroUnit{}, false60 }61 for _, unit := range h.units {62 if pointInRings(lon, lat, unit.rings) {63 return unit, true64 }65 }66 return hydroUnit{}, false67}6869// pointInRings is even-odd ray casting over MultiPolygon nesting: crossing an70// outer ring enters ground, crossing a hole leaves it again.71func pointInRings(lon, lat float64, polygons [][][][]float64) bool {72 inside := false73 for _, polygon := range polygons {74 for _, ring := range polygon {75 for at := 0; at+1 < len(ring); at++ {76 a, b := ring[at], ring[at+1]77 if len(a) < 2 || len(b) < 2 {78 continue79 }80 if (a[1] > lat) == (b[1] > lat) {81 continue82 }83 crossing := a[0] + (lat-a[1])/(b[1]-a[1])*(b[0]-a[0])84 if lon < crossing {85 inside = !inside86 }87 }88 }89 }90 return inside91}9293// featureAnchor is one row's representative position in true degrees: the94// centroid of its largest outer ring, or the middle vertex of its longest line.95// A point row is its own anchor, though no point dataset zones.96func featureAnchor(g geometry) (lon, lat float64, ok bool) {97 switch g.Type {98 case geometryPoint:99 if len(g.Point) < 2 {100 return 0, 0, false101 }102 return g.Point[0], g.Point[1], true103 case geometryRings:104 best := 0.0105 for _, polygon := range g.Rings {106 if len(polygon) == 0 {107 continue108 }109 x, y, area, held := ringCentroid(polygon[0])110 if held && area > best {111 best, lon, lat, ok = area, x, y, true112 }113 }114 return lon, lat, ok115 case geometryLines:116 var longest []float64117 best := 0.0118 for _, line := range g.Lines {119 length := 0.0120 for at := 0; at+1 < len(line); at++ {121 if len(line[at]) < 2 || len(line[at+1]) < 2 {122 continue123 }124 length += math.Hypot(line[at+1][0]-line[at][0], line[at+1][1]-line[at][1])125 }126 if length > best && len(line) > 0 {127 best, longest = length, line[len(line)/2]128 }129 }130 if len(longest) < 2 {131 return 0, 0, false132 }133 return longest[0], longest[1], true134 }135 return 0, 0, false136}137138// ringCentroid is the shoelace centroid of one ring, with the enclosed area to139// rank rings by.140func ringCentroid(ring [][]float64) (lon, lat, area float64, ok bool) {141 var doubled, sumX, sumY float64142 for at := 0; at+1 < len(ring); at++ {143 a, b := ring[at], ring[at+1]144 if len(a) < 2 || len(b) < 2 {145 return 0, 0, 0, false146 }147 cross := a[0]*b[1] - b[0]*a[1]148 doubled += cross149 sumX += (a[0] + b[0]) * cross150 sumY += (a[1] + b[1]) * cross151 }152 if doubled == 0 {153 return 0, 0, 0, false154 }155 return sumX / (3 * doubled), sumY / (3 * doubled), math.Abs(doubled) / 2, true156}157158// hydroClaims tracks, per zone, the subwatersheds its rows' ground lands in,159// and whether any of it landed nowhere.160type hydroClaims struct {161 codes map[string]bool162 missed bool163}164165func (c *hydroClaims) observe(index *hydroIndex, g geometry) {166 for _, position := range claimPositions(g) {167 unit, found := index.locate(position[0], position[1])168 if !found {169 c.missed = true170 continue171 }172 if c.codes == nil {173 c.codes = map[string]bool{}174 }175 c.codes[unit.code] = true176 }177}178179// claimSamples is how many boundary positions a row volunteers beside its180// anchor. The anchor alone would let a zone made of one citywide polygon claim181// whichever subwatershed its centroid happens to sit in; sampling the boundary182// catches the spread, and the cost stays flat no matter how detailed the ring.183const claimSamples = 16184185// claimPositions is the ground a row answers for: its anchor, and an even186// sample of its outer-ring or line vertices.187func claimPositions(g geometry) [][2]float64 {188 var out [][2]float64189 if lon, lat, ok := featureAnchor(g); ok {190 out = append(out, [2]float64{lon, lat})191 }192 var boundary [][]float64193 switch g.Type {194 case geometryRings:195 for _, polygon := range g.Rings {196 if len(polygon) > 0 {197 boundary = append(boundary, polygon[0]...)198 }199 }200 case geometryLines:201 for _, line := range g.Lines {202 boundary = append(boundary, line...)203 }204 }205 step := max((len(boundary)+claimSamples-1)/claimSamples, 1)206 for at := 0; at < len(boundary); at += step {207 if len(boundary[at]) >= 2 {208 out = append(out, [2]float64{boundary[at][0], boundary[at][1]})209 }210 }211 return out212}213214// apply writes the membership onto a zone when every row agreed on one215// subwatershed: a short sentence for the card, and the code as the216// machine-readable key beside it.217func (c *hydroClaims) apply(index *hydroIndex, zone *doc.Feature) {218 if c.missed || len(c.codes) != 1 {219 return220 }221 var code string222 for held := range c.codes {223 code = held224 }225 for _, unit := range index.units {226 if unit.code != code {227 continue228 }229 zone.Description = scrub(fmt.Sprintf("Lies in the %s subwatershed (HUC %s).", unit.name, unit.code))230 if zone.Attrs == nil {231 zone.Attrs = map[string]string{}232 }233 zone.Attrs[semconv.KeyHydroHUC12] = unit.code234 return235 }236}