OpsLog/cmd/ctzipgen/main.go

package main

// One-shot generator: emits internal/uls/ctcounty_gen.go, the Connecticut
// ZIP → legal county table. Not part of the build.
//
// Why it has to exist
// -------------------
// internal/uls resolves a US callsign to a county through GeoNames' ZIP table.
// The Census replaced Connecticut's eight counties with nine PLANNING REGIONS
// as county-equivalents in 2022, and GeoNames followed: 418 of Connecticut's
// 429 ZIPs now report "Capitol Region", "Lower Connecticut River Valley" and
// so on. CQ's USA-CA award did not follow, callbooks did not follow, and no
// operator's log did either — so every Connecticut station resolved to a name
// nothing could match, showed as a new county for ever, and counted for no
// award. Roughly 15 000 US callsigns.
//
// A planning region is NOT a renamed county — it is built from towns drawn
// from several different counties — so there is no name-to-name mapping to be
// had. The ZIP is the only handle, hence this table.
//
// Sources (both public, both free)
// --------------------------------
//	geonames  US.txt from https://download.geonames.org/export/zip/US.zip
//	          — the ZIP list itself, and each ZIP's centroid.
//	crosswalk https://www2.census.gov/geo/docs/maps-data/data/rel2020/
//	          zcta520/tab20_zcta520_county20_natl.txt
//	          — the 2020 ZCTA↔county relationship file, which predates the
//	            change and therefore still carries the eight real counties.
//
// The crosswalk covers 278 of the 418 affected ZIPs. The rest are PO-box-only
// ZIPs, which have no ZCTA at all; each is given the county of the nearest
// resolved ZIP centroid. That is sound here because a PO-box ZIP sits inside
// the town it serves, and Connecticut's counties are tens of kilometres across
// — the fallback can only err on a ZIP that straddles a county line, which the
// ZIP-to-county approach is already documented as accepting (~98%).
//
// Usage:
//
//	go run ./cmd/ctzipgen US.txt tab20_zcta520_county20_natl.txt > internal/uls/ctcounty_gen.go

import (
	"bufio"
	"fmt"
	"math"
	"os"
	"sort"
	"strconv"
	"strings"
)

type zipPt struct {
	zip      string
	lat, lon float64
	county   string // from the crosswalk, "" if the ZIP has no ZCTA
	geoName  string // what GeoNames says today
}

func main() {
	if len(os.Args) < 3 {
		fmt.Fprintln(os.Stderr, "usage: ctzipgen US.txt tab20_zcta520_county20_natl.txt")
		os.Exit(2)
	}

	// 1. Every Connecticut ZIP, with its centroid, from GeoNames.
	var pts []*zipPt
	byZip := map[string]*zipPt{}
	f, err := os.Open(os.Args[1])
	must(err)
	sc := bufio.NewScanner(f)
	sc.Buffer(make([]byte, 0, 64*1024), 256*1024)
	for sc.Scan() {
		c := strings.Split(sc.Text(), "\t")
		if len(c) < 11 || strings.ToUpper(strings.TrimSpace(c[4])) != "CT" {
			continue
		}
		z := strings.TrimSpace(c[1])
		if z == "" || byZip[z] != nil {
			continue
		}
		p := &zipPt{zip: z, lat: atof(c[9]), lon: atof(c[10]), geoName: strings.TrimSpace(c[5])}
		byZip[z] = p
		pts = append(pts, p)
	}
	must(sc.Err())
	f.Close()

	// 2. The 2020 county for each ZCTA. A ZCTA can straddle a county line, so
	//    keep the county holding the largest share of its land area.
	best := map[string]float64{}
	f2, err := os.Open(os.Args[2])
	must(err)
	sc2 := bufio.NewScanner(f2)
	sc2.Buffer(make([]byte, 0, 64*1024), 1<<20)
	sc2.Scan() // header
	for sc2.Scan() {
		c := strings.Split(sc2.Text(), "|")
		if len(c) < 18 || !strings.HasPrefix(c[9], "09") {
			continue
		}
		p := byZip[strings.TrimSpace(c[1])]
		if p == nil {
			continue
		}
		area := atof(c[16])
		if p.county != "" && area <= best[p.zip] {
			continue
		}
		best[p.zip] = area
		p.county = strings.TrimSuffix(strings.TrimSpace(c[10]), " County")
	}
	must(sc2.Err())
	f2.Close()

	// 3. Fill the PO-box-only ZIPs from the nearest ZIP the crosswalk resolved.
	var anchors []*zipPt
	for _, p := range pts {
		if p.county != "" {
			anchors = append(anchors, p)
		}
	}
	if len(anchors) == 0 {
		fmt.Fprintln(os.Stderr, "ctzipgen: crosswalk resolved nothing — wrong file?")
		os.Exit(1)
	}
	filled := 0
	for _, p := range pts {
		if p.county != "" {
			continue
		}
		nearest, bestD := "", math.MaxFloat64
		for _, a := range anchors {
			if d := haversine(p.lat, p.lon, a.lat, a.lon); d < bestD {
				bestD, nearest = d, a.county
			}
		}
		p.county = nearest
		filled++
	}

	// 4. Self-check: the ZIPs GeoNames still labels with a real county must
	//    agree with what we derived, or the derivation is wrong.
	real := map[string]bool{
		"Fairfield": true, "Hartford": true, "Litchfield": true, "Middlesex": true,
		"New Haven": true, "New London": true, "Tolland": true, "Windham": true,
	}
	checked, bad := 0, 0
	for _, p := range pts {
		if !real[p.geoName] {
			continue
		}
		checked++
		if p.geoName != p.county {
			bad++
			fmt.Fprintf(os.Stderr, "MISMATCH %s: geonames %q, derived %q\n", p.zip, p.geoName, p.county)
		}
	}
	fmt.Fprintf(os.Stderr, "ctzipgen: %d zips, %d from crosswalk, %d by nearest; cross-check %d/%d agree\n",
		len(pts), len(pts)-filled, filled, checked-bad, checked)
	if bad > 0 {
		os.Exit(1)
	}

	sort.Slice(pts, func(i, j int) bool { return pts[i].zip < pts[j].zip })

	var b strings.Builder
	b.WriteString("// Code generated by cmd/ctzipgen. DO NOT EDIT.\n\n")
	b.WriteString("package uls\n\n")
	b.WriteString("// ctCounty maps a Connecticut ZIP to its legal county.\n")
	b.WriteString("//\n")
	b.WriteString("// GeoNames reports Connecticut's 2022 planning regions instead, which no\n")
	b.WriteString("// award, callbook or log uses. See cmd/ctzipgen for how this was built and\n")
	b.WriteString("// why a name-to-name mapping cannot work.\n")
	fmt.Fprintf(&b, "var ctCounty = map[string]string{\n")
	for _, p := range pts {
		fmt.Fprintf(&b, "\t%q: %q,\n", p.zip, p.county)
	}
	b.WriteString("}\n")
	fmt.Print(b.String())
}

func atof(s string) float64 {
	v, _ := strconv.ParseFloat(strings.TrimSpace(s), 64)
	return v
}

// haversine returns the great-circle distance in km. Connecticut is small
// enough that a flat approximation would do, but this costs nothing and cannot
// be wrong near the state's edges.
func haversine(lat1, lon1, lat2, lon2 float64) float64 {
	const r = 6371.0
	dLat := (lat2 - lat1) * math.Pi / 180
	dLon := (lon2 - lon1) * math.Pi / 180
	a := math.Sin(dLat/2)*math.Sin(dLat/2) +
		math.Cos(lat1*math.Pi/180)*math.Cos(lat2*math.Pi/180)*math.Sin(dLon/2)*math.Sin(dLon/2)
	return 2 * r * math.Asin(math.Sqrt(a))
}

func must(err error) {
	if err != nil {
		fmt.Fprintln(os.Stderr, "ctzipgen:", err)
		os.Exit(1)
	}
}