1package arcgismap23import (4 "fmt"5 "math"67 "github.com/FelineStateMachine/atlas/internal/mgdoc"8 "github.com/FelineStateMachine/atlas/internal/semconv"9)1011// 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.2021// hydroUnit is one subwatershed as captured: its code, its name, and its22// ground in true degrees.23type hydroUnit struct {24 code string25 name string26 rings [][][][]float6427}2829// hydroIndex answers which subwatershed a point lies in.30type hydroIndex struct {31 units []hydroUnit32}3334// buildHydroIndex reads the captured subwatersheds, if the capture carries35// them; nil means the capture predates the national enrichment and no zone36// makes any claim.37func 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}5657// 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.60func (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}7172// pointInRings is even-odd ray casting over MultiPolygon nesting: crossing73// an outer ring enters ground, crossing a hole leaves it again.74func 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}9596// 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.99func 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}140141// ringCentroid is the shoelace centroid of one ring, with the enclosed142// area to rank rings by.143func 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}160161// hydroClaims tracks, per zone, the subwatersheds its features' ground162// lands in, and whether any of it landed nowhere.163type hydroClaims struct {164 codes map[string]bool165 missed bool166}167168func (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}181182// 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.187const claimSamples = 16188189// claimPositions is the ground a feature answers for: its anchor, and an190// even sample of its outer-ring or line vertices.191func 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}217218// 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.221func (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}