1// Package national joins a volume's own ground to the national hydrography it2// was captured beside.3//4// A city's zones and trails are curated by the city. The subwatersheds under5// them are surveyed by the USGS, and a capture that carries both can say, of6// each piece of the city's ground, which subwatershed it lies in. That sentence7// is worth one line on a card and one machine-readable key beside it, and it is8// computed here rather than at capture time so that re-curating the sentence9// never needs a re-fetch: the evidence base travels in the capture, and the join10// re-runs against it.11//12// # The under-claiming ethic13//14// The claim is made only when every piece of ground a feature answers for15// agrees on one subwatershed. A trail that crosses a divide says nothing. A16// zone whose boundary leaves the surveyed extent says nothing. This is the17// lesson an excised zoning-rules join left behind: one earned sentence beats a18// wall of uncorrelated prose, and a claim that is merely plausible is worse than19// silence because a reader cannot tell the two apart.20package national2122import (23 "encoding/json"24 "fmt"25 "log/slog"26 "math"2728 "github.com/FelineStateMachine/atlas/format/semconv"29 "github.com/FelineStateMachine/atlas/internal/enrich"30 "github.com/FelineStateMachine/atlas/internal/logging"31)3233// Name is what curation queues this enricher as.34const Name = "national"3536// EvidenceName is what the evidence base calls the surveyed subwatersheds.37const EvidenceName = "hydro/huc12.json"3839// The evidence document's identity. Its schema is the enrich lane's, not a40// source's: whoever captured the survey writes this, and this package reads it,41// so the join re-runs from the archive without either side knowing the other's42// vocabulary.43const (44 EvidenceDoc = "atlas-enrich-evidence"45 EvidenceVersion = 146 // EvidenceKind names what this particular evidence is of.47 EvidenceKind = "hydro.huc12"48 // SpaceWorld is the one coordinate space the evidence may be published in:49 // the volume's own projection, which is the space a claim is decided in.50 // Whoever writes the evidence has already projected the survey into the51 // world it was clipped to, so nothing here reprojects anything.52 SpaceWorld = "world"53)5455// Evidence is the surveyed hydrologic units, as the capture carries them.56type Evidence struct {57 Evidence string `json:"evidence"`58 Version int `json:"version"`59 Kind string `json:"kind"`60 Space string `json:"space"`61 Units []Unit `json:"units"`62}6364// Unit is one hydrologic unit: its twelve-digit code, its name, and its ground65// as polygon rings in the volume's own projection.66type Unit struct {67 Code string `json:"code"`68 Name string `json:"name"`69 Rings [][][][2]float64 `json:"rings"`70}7172// ClaimSamples is how many boundary positions a feature volunteers beside its73// anchor, per piece of geometry. The anchor alone would let a zone made of one74// citywide polygon claim whichever subwatershed its centroid happens to sit in;75// sampling the boundary catches the spread, and the cost stays flat no matter76// how detailed the ring.77const ClaimSamples = 167879// Enricher writes the membership join.80type Enricher struct{}8182// New builds the enricher.83func New() *Enricher { return &Enricher{} }8485func (*Enricher) Name() string { return Name }8687func (*Enricher) Declares() []string { return []string{semconv.KeyHydroHUC12} }8889// Enrich reads the evidence base and claims what it can be sure of.90//91// An evidence base that carries nothing about this volume is not an error: it92// is a volume nobody surveyed, and the enricher says nothing about it.93func (e *Enricher) Enrich(v *enrich.Volume, ctx enrich.Context) (enrich.Contribution, error) {94 out := enrich.Contribution{Enricher: Name, Volume: v.Slug}95 if ctx.Evidence == nil {96 return out, nil97 }98 log := ctx.Log99 if log == nil {100 log = slog.New(slog.DiscardHandler)101 }102 log = log.With(logging.Op(Name), logging.Volume(v.Slug))103104 data, held, err := ctx.Evidence.Open(EvidenceName)105 if err != nil {106 return enrich.Contribution{}, fmt.Errorf("read %s: %w", EvidenceName, err)107 }108 if !held {109 log.Debug("no hydrography was captured beside this volume")110 return out, nil111 }112 evidence, err := ReadEvidence(data)113 if err != nil {114 return enrich.Contribution{}, err115 }116 if len(evidence.Units) == 0 {117 return out, nil118 }119120 // The evidence does not make claims about itself: the hydrography is the121 // survey, and a subwatershed needs no sentence saying which subwatershed it122 // lies in.123 surveyed := make(map[string]bool)124 for _, title := range ctx.Curation.EvidenceCollections(v.Slug) {125 surveyed[title] = true126 }127128 claimed, straddled := 0, 0129 for _, world := range v.Worlds {130 for _, collection := range world.Collections {131 if collection.Kind == enrich.KindPoint || surveyed[collection.Title] {132 continue133 }134 for _, feature := range collection.Features {135 if feature.Attrs[semconv.KeyHydroHUC12] != "" {136 continue137 }138 unit, sure := claim(evidence.Units, feature)139 if !sure {140 straddled++141 continue142 }143 claimed++144 out.Ops = append(out.Ops, enrich.Op{145 Kind: enrich.OpSetAttr,146 World: world.Slug,147 Feature: feature.ID,148 Entity: semconv.EntityFeature,149 Key: semconv.KeyHydroHUC12,150 Value: unit.Code,151 })152 if feature.Description == "" {153 out.Ops = append(out.Ops, enrich.Op{154 Kind: enrich.OpSetProse,155 World: world.Slug,156 Feature: feature.ID,157 Value: Sentence(unit),158 })159 }160 }161 }162 }163 if claimed > 0 || straddled > 0 {164 log.Info("hydrologic membership joined", "claimed", claimed, "unclaimed", straddled,165 "units", len(evidence.Units))166 }167 return out, nil168}169170// Sentence is the one line a claim earns on a card.171func Sentence(unit Unit) string {172 return fmt.Sprintf("Lies in the %s subwatershed (HUC %s).", unit.Name, unit.Code)173}174175// ReadEvidence decodes an evidence document and holds it to what this package176// can actually reason about.177func ReadEvidence(data []byte) (Evidence, error) {178 var evidence Evidence179 if err := json.Unmarshal(data, &evidence); err != nil {180 return Evidence{}, fmt.Errorf("decode %s: %w", EvidenceName, err)181 }182 if evidence.Evidence != EvidenceDoc {183 return Evidence{}, fmt.Errorf("evidence says %q, not %q", evidence.Evidence, EvidenceDoc)184 }185 if evidence.Version != EvidenceVersion {186 return Evidence{}, fmt.Errorf("evidence version %d, want %d", evidence.Version, EvidenceVersion)187 }188 if evidence.Kind != EvidenceKind {189 return Evidence{}, fmt.Errorf("evidence is of %q, not %q", evidence.Kind, EvidenceKind)190 }191 if evidence.Space != SpaceWorld {192 return Evidence{}, fmt.Errorf(193 "evidence is published in %q; this join reads %q, the volume's own projection",194 evidence.Space, SpaceWorld)195 }196 for _, unit := range evidence.Units {197 if err := semconv.Check(semconv.EntityFeature, semconv.KeyHydroHUC12, unit.Code); err != nil {198 return Evidence{}, fmt.Errorf("evidence unit %q: %w", unit.Name, err)199 }200 }201 return evidence, nil202}203204// claim answers which unit a feature's ground lies in, and whether the answer205// is unanimous. Every position the feature answers for has to land in the same206// unit, and none may land outside the survey.207func claim(units []Unit, feature enrich.Feature) (Unit, bool) {208 found := -1209 for _, geometry := range feature.Geometry {210 for _, position := range Positions(geometry) {211 at := locate(units, position[0], position[1])212 if at < 0 {213 return Unit{}, false214 }215 if found >= 0 && at != found {216 return Unit{}, false217 }218 found = at219 }220 }221 if found < 0 {222 return Unit{}, false223 }224 return units[found], true225}226227// locate answers which unit holds a position, by index, or -1. Units are228// visited in the order the evidence lists them, so a position exactly on a229// shared boundary answers the same way every run.230func locate(units []Unit, lng, lat float64) int {231 for index, unit := range units {232 if inRings(unit.Rings, lng, lat) {233 return index234 }235 }236 return -1237}238239// inRings is even-odd ray casting over polygon nesting: crossing an outer ring240// enters ground, crossing a hole leaves it again.241func inRings(polygons [][][][2]float64, lng, lat float64) bool {242 inside := false243 for _, polygon := range polygons {244 for _, ring := range polygon {245 for at := 0; at+1 < len(ring); at++ {246 a, b := ring[at], ring[at+1]247 if (a[1] > lat) == (b[1] > lat) {248 continue249 }250 crossing := a[0] + (lat-a[1])/(b[1]-a[1])*(b[0]-a[0])251 if lng < crossing {252 inside = !inside253 }254 }255 }256 }257 return inside258}259260// Positions is the ground one piece of a feature's geometry answers for: its261// anchor, and an even sample of its outer-ring or line vertices.262//263// It is exported because it is the whole of what "how much of a feature is264// looked at" means, and a test that wants to say what a claim stood on has to265// be able to ask.266func Positions(g enrich.Geometry) [][2]float64 {267 var out [][2]float64268 var boundary [][2]float64269 best := 0.0270 var anchor [2]float64271 held := false272273 for _, polygon := range g.Rings() {274 if len(polygon) == 0 {275 continue276 }277 if lng, lat, area, ok := ringCentroid(polygon[0]); ok && area > best {278 best, anchor, held = area, [2]float64{lng, lat}, true279 }280 boundary = append(boundary, polygon[0]...)281 }282 for _, line := range g.Lines() {283 length := 0.0284 for at := 0; at+1 < len(line); at++ {285 length += math.Hypot(line[at+1][0]-line[at][0], line[at+1][1]-line[at][1])286 }287 if length > best && len(line) > 0 {288 best, anchor, held = length, line[len(line)/2], true289 }290 boundary = append(boundary, line...)291 }292 if !held {293 // A point feature is its own anchor, though no point collection is294 // joined today.295 positions := g.Positions()296 if len(positions) == 1 {297 return positions298 }299 return nil300 }301 out = append(out, anchor)302 step := max((len(boundary)+ClaimSamples-1)/ClaimSamples, 1)303 for at := 0; at < len(boundary); at += step {304 out = append(out, boundary[at])305 }306 return out307}308309// ringCentroid is the shoelace centroid of one ring, with the enclosed area to310// rank rings by.311func ringCentroid(ring [][2]float64) (lng, lat, area float64, ok bool) {312 var doubled, sumX, sumY float64313 for at := 0; at+1 < len(ring); at++ {314 a, b := ring[at], ring[at+1]315 cross := a[0]*b[1] - b[0]*a[1]316 doubled += cross317 sumX += (a[0] + b[0]) * cross318 sumY += (a[1] + b[1]) * cross319 }320 if doubled == 0 {321 return 0, 0, 0, false322 }323 return sumX / (3 * doubled), sumY / (3 * doubled), math.Abs(doubled) / 2, true324}