Compare commits

...
4 Commits
Author SHA1 Message Date
Knou a9e722791d Complete Köppen-Geiger filter review with Mixed climate class 2026-09-15 14:32:15 -04:00
KnouandClaude Opus 5 4d2b3e3d44 Complete Köppen-Geiger filter review with Mixed climate class
Classify each county by area-weighted Köppen class shares: a county is
predominantly its top class when that class covers at least 50% of its
land and leads the runner-up by at least 5 percentage points; otherwise
it is Mixed (133 of 3,143 counties in the 50 states and DC).

- Add build_county_koppen_metric.py (writes data/metrics/koppen.csv) and
  apply_koppen_metric_to_climate_data.py (writes koppenZone plus
  koppenPrimaryClass/koppenSecondaryClass for Mixed counties).
- Move shared helpers into scripts/common/ (county loading, Köppen
  legend, area-weighted raster shares); fix the 180th-meridian raster
  window for Aleutians West.
- Add check_climate_data.py to validate the app CSV.
- Draw Mixed counties in app.js as diagonal stripes of their top two
  classes, fixed to the ground and following the map at every zoom, with
  a crossfade only when the stripe size changes. Filtering a class also
  matches Mixed counties where it is primary or secondary.
- Document the rule, display, and pipeline plan in docs/ and update the
  README and data-source notes.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
2026-09-14 02:54:25 -04:00
Knou 92fbbfb2e9 README.md has been updated to reflect project changes. 2026-09-11 16:30:15 -04:00
Knou 4e2d878e5b Refine climate metrics and data pipeline 2026-07-24 22:55:14 -04:00
47 changed files with 6867 additions and 3650 deletions
+1
View File
@@ -22,3 +22,4 @@ data/*
# Runtime logs # Runtime logs
*.log *.log
/.mcp.json
+86 -14
View File
@@ -14,7 +14,7 @@ the climate side of the analysis.
- Interactive Leaflet map with county boundaries and county selection - Interactive Leaflet map with county boundaries and county selection
- Metric groups for climate classification, temperature, precipitation, - Metric groups for climate classification, temperature, precipitation,
moisture, and solar exposure moisture, and solar resource
- Numeric range filters and categorical filters - Numeric range filters and categorical filters
- County detail panel showing all available climate metrics - County detail panel showing all available climate metrics
- Legends and in-app source descriptions - Legends and in-app source descriptions
@@ -28,15 +28,18 @@ The browser currently loads 3,221 county-level records from
| Group | Metrics | | Group | Metrics |
| --- | --- | | --- | --- |
| Koppen-Geiger Classification | Majority county climate class | | Koppen-Geiger Classification | Predominant county climate class (covering at least 50% of the county's land and leading the runner-up by at least 5 points), or Mixed, drawn as stripes of the county's top two classes |
| Temperature & Extremes | Annual average temperature, diurnal temperature range, absolute extreme days per year, and heat-index days | | Temperature & Extremes | Annual average temperature, diurnal temperature range, annual extreme temperature days, and annual 90 F+ heat-index days |
| Precipitation & Moisture | Annual precipitation, precipitation seasonality, wettest month, driest month, and summer specific humidity | | Precipitation & Moisture | Annual precipitation, precipitation seasonality, wettest month, driest month, and summer specific humidity |
| Solar Exposure | Average daily global horizontal irradiance (GHI) and cloudiness index | | Solar Resource | Mean daily global horizontal radiation (GHI) and clear-sky GHI reduction index |
Most long-term climate metrics use a 1991-2020 reference period. The absolute Most long-term climate metrics use a 1991-2020 reference period. In the
extreme-days metric currently uses 1991-2025 data. Definitions and time periods checked-in CSV, the annual extreme-temperature-days metric uses 1991-2025 data
are shown in the application's **Sources** dialog and documented in more detail and is the average annual count of days with `Tmax >= 95 F` or `Tmin <= 0 F`.
in [`scripts/county_data_sources.md`](scripts/county_data_sources.md). The build script defaults its analysis end year to the latest likely complete
calendar year. Definitions and time periods are shown in the application's
**Sources** dialog and documented in more detail in
[`scripts/county_data_sources.md`](scripts/county_data_sources.md).
## Run the Explorer ## Run the Explorer
@@ -89,9 +92,13 @@ will not work.
| `app.js` | Map rendering, filtering, county details, and source metadata | | `app.js` | Map rendering, filtering, county details, and source metadata |
| `serve.ps1` | Local HTTP server launcher | | `serve.ps1` | Local HTTP server launcher |
| `data/climate-data.csv` | Browser-ready county climate records | | `data/climate-data.csv` | Browser-ready county climate records |
| `data/metrics/` | Per-metric county outputs, starting with `koppen.csv` |
| `data/geojson-counties-fips.json` | County geometry keyed by FIPS code | | `data/geojson-counties-fips.json` | County geometry keyed by FIPS code |
| `scripts/` | Climate-data download, aggregation, and update tools | | `scripts/` | Climate-data download, aggregation, and update tools |
| `tests/` | Tests for the NSRDB polygon request and download workflow | | `scripts/county_data_sources.md` | Detailed metric definitions, data provenance, and pipeline examples |
| `scripts/requirements_county_etl.txt` | Python dependencies for the offline data pipeline |
| `docs/` | Filter calculations, the pipeline plan, and design notes |
| `tests/` | Tests for the Köppen metric, the climate CSV check, NSRDB request/download wrappers, cloud-metric merging, and gridMET heat-index calculations |
## Climate Data Pipeline ## Climate Data Pipeline
@@ -103,10 +110,10 @@ outputs from sources including:
- NOAA NCEI Climate Normals and nClimGrid data - NOAA NCEI Climate Normals and nClimGrid data
- gridMET humidity data - gridMET humidity data
- NREL National Solar Radiation Database data - NREL National Solar Radiation Database data
- U.S. Census Bureau county geometry - Plotly's county GeoJSON, keyed by U.S. Census county FIPS codes
To work on the Python data pipeline, create a virtual environment and install The Python tooling is configured for Python 3.11. To work on the data pipeline,
the ETL dependencies: create a virtual environment and install the ETL dependencies:
```powershell ```powershell
python -m venv .venv python -m venv .venv
@@ -118,6 +125,71 @@ Some data-generation workflows download large files or require NREL/NSRDB API
credentials. Generated source datasets are intentionally excluded from Git; credentials. Generated source datasets are intentionally excluded from Git;
only the browser-ready CSV and GeoJSON are tracked. only the browser-ready CSV and GeoJSON are tracked.
Run script commands from the project root so their default `data/...` paths
resolve correctly. The pipeline is split into a base generator, source-specific
builders, and small scripts that merge the resulting metrics into
`data/climate-data.csv`. Most `apply_*.py` commands update that CSV in place by
default.
### Script Inventory
| Script | Current role |
| --- | --- |
| `build_county_climate_data.py` | Builds the base app CSV from county geometry, Koppen-Geiger data, NOAA temperature/precipitation data, and optional solar inputs. Its older `extremeDays` output is replaced by the current absolute-threshold stage below, and its largest-share `koppenZone` by the Köppen stage. It writes only the base columns, so do not run it over the live CSV. |
| `build_county_koppen_metric.py` / `apply_koppen_metric_to_climate_data.py` | Builds area-weighted Köppen class shares per county in `data/metrics/koppen.csv`, classifies each county as predominant or Mixed, and writes `koppenZone` plus the two stripe-class columns for Mixed counties. |
| `apply_precipitation_month_metrics_to_climate_data.py` | Recomputes and merges the 1991-2020 wettest- and driest-month categories from monthly nClimGrid precipitation. |
| `build_county_locally_extreme_data.py` | Downloads or reads cached nClimGrid-Daily county Tmax/Tmin files, calculates county-percentile diagnostics, and calculates the app-facing absolute 95 F / 0 F day counts. |
| `apply_locally_extreme_metric_to_climate_data.py` | Writes `absoluteExtremeDays`, removes retired locally extreme/legacy fields, and selects polygon GHI with representative-point GHI as fallback. The filename is retained from the earlier pipeline. |
| `build_county_diurnal_temperature_range.py` / `apply_diurnal_temperature_range_to_climate_data.py` | Builds the 1991-2020 county mean daily Tmax-minus-Tmin artifact and merges `avgDiurnalTempRangeF`. |
| `download_gridmet_data.py` | Downloads 1991-2020 `sph`, `rmax`, and `rmin` NetCDF files by default. |
| `summarize_county_gridmet_humidity.py` / `apply_gridmet_humidity_metric_to_climate_data.py` | Produces county summer specific humidity and a 90 F+ Heat Index day proxy, then merges those metrics and FIPS audit fields. |
| `build_county_representative_points.py` | Creates interior county points used by the lightweight NSRDB workflows. |
| `fetch_nsrdb_representative_point_ghi.py` / `rebuild_nsrdb_representative_point_ghi_summary.py` | Fetches point-based NSRDB GHI or rebuilds its summary from cached responses without another API call. |
| `fetch_nsrdb_representative_point_cloud_metrics.py` | Fetches point-based GHI, clear-sky GHI, and cloud type, then summarizes the clear-sky GHI reduction index. |
| `request_nsrdb_county_polygon_archives.py` | Shared, resumable NSRDB polygon-request engine with tiling, site-count checks, pacing, and Polar fallback. |
| `request_nsrdb_county_polygon_ghi_archives.py` / `request_nsrdb_county_polygon_cloud_archives.py` | Recommended wrappers around the shared request engine, with separate GHI and cloud attributes and output paths. |
| `download_nsrdb_county_polygon_archives.py` | Shared state-machine downloader for completed NSRDB archive jobs. |
| `download_nsrdb_county_polygon_ghi_archives.py` / `download_nsrdb_county_polygon_cloud_archives.py` | Recommended wrappers around the shared downloader, keeping GHI and cloud archives separate. |
| `summarize_nsrdb_county_polygon_archives.py` | Combines county/tile GHI archives into area-weighted county summaries. |
| `summarize_nsrdb_county_polygon_cloud_archives.py` | Combines county/tile cloud archives into area-weighted clear-sky GHI reduction summaries. |
| `apply_nsrdb_cloud_metric_to_climate_data.py` | Merges the clear-sky GHI reduction index, preferring polygon summaries and falling back to representative points. |
| `check_climate_data.py` | Validates `data/climate-data.csv`, including Köppen codes and the stripe-class columns. |
| `common/` | Shared helpers: county loading, the Köppen legend, and area-weighted raster shares. |
The metric-specific NSRDB request and download wrappers are the normal entry
points. The shared engines remain available for custom attributes or artifact
paths. Representative-point results provide a faster first pass; polygon
summaries are the preferred county-area result when available.
Common local enrichment stages, after their source files have been downloaded,
are:
```powershell
.venv\Scripts\python.exe scripts\build_county_koppen_metric.py
.venv\Scripts\python.exe scripts\apply_koppen_metric_to_climate_data.py
.venv\Scripts\python.exe scripts\build_county_locally_extreme_data.py --skip-download
.venv\Scripts\python.exe scripts\apply_locally_extreme_metric_to_climate_data.py
.venv\Scripts\python.exe scripts\build_county_diurnal_temperature_range.py
.venv\Scripts\python.exe scripts\apply_diurnal_temperature_range_to_climate_data.py
.venv\Scripts\python.exe scripts\summarize_county_gridmet_humidity.py --years 1991-2020
.venv\Scripts\python.exe scripts\apply_gridmet_humidity_metric_to_climate_data.py
.venv\Scripts\python.exe scripts\apply_nsrdb_cloud_metric_to_climate_data.py
```
The order matters when rebuilding from scratch: run the Köppen apply step after
the base build, generate the locally extreme comparison and solar summaries
before running their apply step, and summarize gridMET or NSRDB downloads
before merging them. See the data-source document linked above for acquisition
commands, expected artifacts, FIPS handling, and the representative-point and
polygon NSRDB workflows. The full rebuild order is in
[`docs/pipeline-plan.md`](docs/pipeline-plan.md).
After updating the CSV, validate it with:
```powershell
.venv\Scripts\python.exe scripts\check_climate_data.py
```
Run the current automated tests with: Run the current automated tests with:
```powershell ```powershell
@@ -128,8 +200,8 @@ python -m unittest discover -s tests
The longer-term goal is to add mood-based metrics and investigate whether The longer-term goal is to add mood-based metrics and investigate whether
patterns in those metrics are associated with climate characteristics such as patterns in those metrics are associated with climate characteristics such as
temperature, sunlight, cloudiness, humidity, precipitation, or extreme-weather temperature, sunlight availability, clear-sky GHI reduction, humidity,
frequency. precipitation, or extreme-weather frequency.
That phase still requires decisions about: That phase still requires decisions about:
+412 -39
View File
@@ -1,6 +1,6 @@
"use strict"; "use strict";
const APP_ASSET_VERSION = "precip-months-cloudiness-20260601"; const APP_ASSET_VERSION = "koppen-mixed-cleanup-20260914";
const US_COUNTIES_GEOJSON_URL = `data/geojson-counties-fips.json?v=${APP_ASSET_VERSION}`; const US_COUNTIES_GEOJSON_URL = `data/geojson-counties-fips.json?v=${APP_ASSET_VERSION}`;
const CLIMATE_DATA_CSV_URL = `data/climate-data.csv?v=${APP_ASSET_VERSION}`; const CLIMATE_DATA_CSV_URL = `data/climate-data.csv?v=${APP_ASSET_VERSION}`;
const NO_DATA_FILL_COLOR = "#d9dde8"; const NO_DATA_FILL_COLOR = "#d9dde8";
@@ -102,9 +102,34 @@ const KOPPEN_CLASS_META = {
Dfc: { label: "Subarctic", color: "#007D7D" }, Dfc: { label: "Subarctic", color: "#007D7D" },
Dfd: { label: "Extremely Cold Subarctic", color: "#00465F" }, Dfd: { label: "Extremely Cold Subarctic", color: "#00465F" },
ET: { label: "Tundra", color: "#B2B2B2" }, ET: { label: "Tundra", color: "#B2B2B2" },
EF: { label: "Ice Cap", color: "#666666" } EF: { label: "Ice Cap", color: "#666666" },
Mixed: { label: "Mixed Climate", color: "#9aa3b2" }
}; };
// Counties with no predominant Koppen-Geiger class; see docs/koppen-mixed-display-plan.md.
const KOPPEN_MIXED_CLASS = "Mixed";
// Standard style: share of each stripe pair given to the primary (top) class; the secondary class
// gets the rest.
const KOPPEN_MIXED_PRIMARY_STRIPE_FRACTION = 0.8;
// Stripe pairs per 100 miles, measured across the stripes at DEFAULT_COUNTRY_VIEW. Stripes then
// follow the map, so a county keeps the same number of stripes at every zoom.
const KOPPEN_MIXED_LINE_PAIRS_PER_100_MILES = 5;
// Below this width (zoomed far out), stripes switch to the far-out style.
const KOPPEN_MIXED_MIN_PAIR_WIDTH_PX = 5;
// Far-out style: the map-following width doubled until it reaches the minimum (fewer lines),
// with a thicker secondary.
const KOPPEN_MIXED_FAR_PRIMARY_STRIPE_FRACTION = 0.7;
// Length of the crossfade when the stripe size changes (the far-out zooms and their switch to standard).
const KOPPEN_MIXED_STYLE_FADE_MS = 250;
// SVG rotation for the stripes; -45 draws them from upper left to lower right.
const KOPPEN_MIXED_STRIPE_ROTATION_DEGREES = -45;
const KOPPEN_MIXED_LEGEND_PRIMARY_GRAY = "#d4d8df";
const KOPPEN_MIXED_LEGEND_SECONDARY_GRAY = "#6b7383";
const METERS_PER_MILE = 1609.344;
// Web Mercator ground resolution at zoom 0 on the equator, for 256-pixel tiles.
const WEB_MERCATOR_METERS_PER_PIXEL_AT_ZOOM_0 = 156543.03392;
const SVG_NAMESPACE = "http://www.w3.org/2000/svg";
const MONTH_CATEGORY_META = { const MONTH_CATEGORY_META = {
January: { label: "January", color: "#355c9a" }, January: { label: "January", color: "#355c9a" },
February: { label: "February", color: "#4d79b8" }, February: { label: "February", color: "#4d79b8" },
@@ -129,7 +154,7 @@ const METRICS = {
allOptionLabel: "All Koppen Classes", allOptionLabel: "All Koppen Classes",
source: { source: {
description: description:
"Beck et al. 2023 Koppen-Geiger maps. County values are generated offline by assigning each county its majority Koppen-Geiger raster class.", "Beck et al. 2023 Koppen-Geiger maps (1991-2020, 1 km). County class shares are area-weighted. A county shows its top class when that class covers at least 50% of its land and leads the runner-up by at least 5 percentage points; otherwise it is Mixed Climate, striped in the colors of its top two classes.",
links: [ links: [
{ {
label: "Beck et al. 2023 Koppen-Geiger", label: "Beck et al. 2023 Koppen-Geiger",
@@ -162,7 +187,6 @@ const METRICS = {
min: 14, min: 14,
max: 33, max: 33,
boundsStep: 0.1, boundsStep: 0.1,
sliderStep: 0.1,
unit: "F", unit: "F",
palette: ["#f2f7f5", "#bcdccf", "#83bba6", "#4b927d", "#1f5f54"], palette: ["#f2f7f5", "#bcdccf", "#83bba6", "#4b927d", "#1f5f54"],
source: { source: {
@@ -247,12 +271,11 @@ const METRICS = {
} }
}, },
absoluteExtremeDays: { absoluteExtremeDays: {
label: "Absolute Extreme Days / Year", label: "Annual Extreme Temperature Days",
type: "numeric", type: "numeric",
min: 0, min: 0,
max: 150, max: 150,
boundsStep: 1, boundsStep: 1,
sliderStep: 1,
unit: "days", unit: "days",
palette: ["#fff9e6", "#ffe49a", "#ffd067", "#e7af20", "#b88409"], palette: ["#fff9e6", "#ffe49a", "#ffd067", "#e7af20", "#b88409"],
source: { source: {
@@ -267,21 +290,24 @@ const METRICS = {
} }
}, },
humidHeatDays: { humidHeatDays: {
label: "Heat Index", label: "Annual 90 F+ Heat Index Days",
type: "numeric", type: "numeric",
min: 0, min: 0,
max: 180, max: 180,
boundsStep: 1, boundsStep: 1,
sliderStep: 0.1,
unit: "days", unit: "days",
palette: ["#fff7ec", "#fdd49e", "#fdbb84", "#ef6548", "#b30000"], palette: ["#fff7ec", "#fdd49e", "#fdbb84", "#ef6548", "#b30000"],
source: { source: {
description: description:
"Heat Index is estimated from NOAA nClimGrid-Daily county Tmax and gridMET daily minimum relative humidity (rmin), 1991-2020. This metric is the average annual count of days where estimated daily Heat Index is at least 90 F; county-vintage joins and proxies are documented in the humidHeatSourceFips and humidHeatFipsAdjustment fields.", "Average annual count of days with an estimated Heat Index of at least 90 F, the lower bound of the NWS Extreme Caution category. The estimate pairs NOAA nClimGrid-Daily county Tmax with gridMET daily minimum relative humidity (rmin), 1991-2020; it is a daily-extrema proxy rather than an observed hourly maximum. County-vintage joins and proxies are documented in the humidHeatSourceFips and humidHeatFipsAdjustment fields.",
links: [ links: [
{ {
label: "NOAA National Weather Service Heat Index", label: "NWS Heat Index Categories",
url: "https://www.weather.gov/safety/heat-index" url: "https://www.weather.gov/ama/heatindex"
},
{
label: "NWS Heat Index Equation",
url: "https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml"
}, },
{ {
label: "NOAA nClimGrid-Daily", label: "NOAA nClimGrid-Daily",
@@ -294,41 +320,39 @@ const METRICS = {
] ]
} }
}, },
avgSolarGhiKwhM2Day: { meanDailyGlobalHorizontalRadiationKwhM2Day: {
label: "Avg Daily Solar Irradiance (GHI)", label: "Mean Daily Global Horizontal Radiation (GHI)",
type: "numeric", type: "numeric",
min: 2.5, min: 2.5,
max: 6.5, max: 6.5,
boundsStep: 0.1, boundsStep: 0.1,
sliderStep: 0.1,
unit: "kWh/m2/day", unit: "kWh/m2/day",
palette: ["#fff8d8", "#f8df83", "#e8b94e", "#c9892b", "#985f1b"], palette: ["#fff8d8", "#f8df83", "#e8b94e", "#c9892b", "#985f1b"],
source: { source: {
description: description:
"NREL National Solar Radiation Database TMY PSM v4. County polygon archive summaries use area-weighted GHI; representative-point values are kept as fallback where polygon summaries are unavailable.", "National Renewable Energy Laboratory (NREL) National Solar Radiation Database (NSRDB) TMY PSM v4. County polygon archive summaries use area-weighted GHI; representative-point values are kept as fallback where polygon summaries are unavailable.",
links: [ links: [
{ {
label: "NREL NSRDB", label: "NREL National Solar Radiation Database",
url: "https://nsrdb.nrel.gov/" url: "https://nsrdb.nrel.gov/"
} }
] ]
} }
}, },
cloudinessIndexPct: { clearSkyGhiReductionIndex: {
label: "Cloudiness Index", label: "Clear-Sky GHI Reduction Index",
type: "numeric", type: "numeric",
min: 0.05, min: 0.05,
max: 0.6, max: 0.6,
boundsStep: 0.01, boundsStep: 0.01,
sliderStep: 0.01,
unit: "ratio", unit: "ratio",
palette: ["#fff7bc", "#d9f0d3", "#a6cee3", "#6baed6", "#4b5563"], palette: ["#fff7bc", "#d9f0d3", "#a6cee3", "#6baed6", "#4b5563"],
source: { source: {
description: description:
"NREL National Solar Radiation Database TMY PSM v4 representative-point summary. The index is 1 minus the daylight average of observed GHI divided by modeled Clearsky GHI, clamped to a 0-1 range; higher values indicate more sunlight reduction relative to clear-sky conditions.", "National Renewable Energy Laboratory (NREL) National Solar Radiation Database (NSRDB) TMY PSM v4. The index is 1 minus the daylight average of observed GHI divided by modeled clear-sky GHI, clamped to a 0-1 range; higher values indicate more solar radiation reduction relative to clear-sky conditions.",
links: [ links: [
{ {
label: "NREL NSRDB", label: "NREL National Solar Radiation Database",
url: "https://nsrdb.nrel.gov/" url: "https://nsrdb.nrel.gov/"
} }
] ]
@@ -340,7 +364,6 @@ const METRICS = {
min: 4, min: 4,
max: 20, max: 20,
boundsStep: 0.1, boundsStep: 0.1,
sliderStep: 0.1,
unit: "g/kg", unit: "g/kg",
palette: ["#f6f0d6", "#d7dfad", "#9cc895", "#58a9a3", "#2e6f8e"], palette: ["#f6f0d6", "#d7dfad", "#9cc895", "#58a9a3", "#2e6f8e"],
source: { source: {
@@ -380,9 +403,9 @@ const METRIC_GROUPS = [
metricKeys: ["annualPrecipIn", "seasonalityIndex", "wettestPrecipMonth", "driestPrecipMonth", "avgSummerSpecificHumidityGKg"] metricKeys: ["annualPrecipIn", "seasonalityIndex", "wettestPrecipMonth", "driestPrecipMonth", "avgSummerSpecificHumidityGKg"]
}, },
{ {
key: "solarExposure", key: "solarResource",
label: "Solar Exposure", label: "Solar Resource",
metricKeys: ["avgSolarGhiKwhM2Day", "cloudinessIndexPct"] metricKeys: ["meanDailyGlobalHorizontalRadiationKwhM2Day", "clearSkyGhiReductionIndex"]
} }
].map((group) => ({ ].map((group) => ({
...group, ...group,
@@ -493,6 +516,9 @@ let didMapDragGesture = false;
let clearTooltipSuppressionOnNextMove = false; let clearTooltipSuppressionOnNextMove = false;
let activeCountyTooltip = null; let activeCountyTooltip = null;
let focusedElementBeforeSourcesModal = null; let focusedElementBeforeSourcesModal = null;
let koppenStripeDefs = null;
let koppenStripeStyle = null;
let koppenStripeFadeFrame = null;
// Find a required DOM element and fail early if the page markup is missing it. // Find a required DOM element and fail early if the page markup is missing it.
function getRequiredElement(id) { function getRequiredElement(id) {
@@ -718,6 +744,10 @@ function normalizeKoppenCode(code) {
return null; return null;
} }
if (trimmed.toLowerCase() === KOPPEN_MIXED_CLASS.toLowerCase()) {
return KOPPEN_MIXED_CLASS;
}
if (!/^[A-Za-z]{2,3}$/.test(trimmed)) { if (!/^[A-Za-z]{2,3}$/.test(trimmed)) {
return null; return null;
} }
@@ -893,6 +923,10 @@ function sanitizeOverrideRecord(rawRecord) {
} }
}); });
// Stripe classes for Mixed Koppen counties; blank for predominant counties.
sanitizedRecord.koppenPrimaryClass = normalizeKoppenCode(rawRecord.koppenPrimaryClass);
sanitizedRecord.koppenSecondaryClass = normalizeKoppenCode(rawRecord.koppenSecondaryClass);
return sanitizedRecord; return sanitizedRecord;
} }
@@ -965,6 +999,17 @@ function roundMetricBound(value, step, direction) {
return Number(rounded.toFixed(decimalPlaces)); return Number(rounded.toFixed(decimalPlaces));
} }
// Choose a power-of-ten slider step that keeps three significant digits of range precision.
function getMetricSliderStep(metric) {
const largestMagnitude = Math.max(Math.abs(metric.min), Math.abs(metric.max));
if (!Number.isFinite(largestMagnitude) || largestMagnitude === 0) {
return 1;
}
const magnitudeExponent = Math.floor(Math.log10(largestMagnitude));
return 10 ** (magnitudeExponent - 2);
}
// Collect the category values that are present in loaded county records for a metric. // Collect the category values that are present in loaded county records for a metric.
function getUniqueCategoryValuesInData(metricKey) { function getUniqueCategoryValuesInData(metricKey) {
const valueSet = new Set(); const valueSet = new Set();
@@ -972,6 +1017,12 @@ function getUniqueCategoryValuesInData(metricKey) {
if (stats[metricKey]) { if (stats[metricKey]) {
valueSet.add(stats[metricKey]); valueSet.add(stats[metricKey]);
} }
// Classes that appear only as stripe colors in Mixed counties are still listed.
const stripeClasses = metricKey === "koppenZone" ? getKoppenStripeClasses(stats) : null;
if (stripeClasses) {
valueSet.add(stripeClasses.primary);
valueSet.add(stripeClasses.secondary);
}
}); });
const categoryOrder = Object.keys(METRICS[metricKey]?.categoryMeta || {}); const categoryOrder = Object.keys(METRICS[metricKey]?.categoryMeta || {});
@@ -1000,7 +1051,7 @@ function getCategoricalDisplayLabel(metricKey, value) {
if (!meta) { if (!meta) {
return value; return value;
} }
return metricKey === "koppenZone" ? `${value} - ${meta.label}` : meta.label; return metricKey === "koppenZone" && value !== KOPPEN_MIXED_CLASS ? `${value} - ${meta.label}` : meta.label;
} }
// Set numeric metric slider bounds from the loaded county data range. // Set numeric metric slider bounds from the loaded county data range.
@@ -1106,6 +1157,298 @@ function colorForValue(metricKey, value) {
return metric.palette[index]; return metric.palette[index];
} }
// Convert the Mixed line pair density into one stripe pair's width in pixels at the default view.
function getKoppenMixedStripePairWidthPx() {
const latitudeRadians = (DEFAULT_COUNTRY_VIEW.center[0] * Math.PI) / 180;
const metersPerPixel =
(WEB_MERCATOR_METERS_PER_PIXEL_AT_ZOOM_0 * Math.cos(latitudeRadians)) / 2 ** DEFAULT_COUNTRY_VIEW.zoom;
const pixelsPer100Miles = (100 * METERS_PER_MILE) / metersPerPixel;
return pixelsPer100Miles / KOPPEN_MIXED_LINE_PAIRS_PER_100_MILES;
}
// Return a Mixed county's primary and secondary Koppen classes, or null when it is not striped.
function getKoppenStripeClasses(stats) {
if (stats?.koppenZone !== KOPPEN_MIXED_CLASS) {
return null;
}
const primary = stats.koppenPrimaryClass;
const secondary = stats.koppenSecondaryClass;
const isStripeClass = (code) => Boolean(KOPPEN_CLASS_META[code]) && code !== KOPPEN_MIXED_CLASS;
return isStripeClass(primary) && isStripeClass(secondary) ? { primary, secondary } : null;
}
// Build the SVG pattern id for one primary/secondary stripe pair.
function getKoppenStripePatternId(primary, secondary) {
return `koppen-mixed-${primary}-${secondary}`;
}
// Return the stripe style for a zoom level. Every width is the reference width times a power of
// two, so the old and new layouts of any zoom change fit one pattern tile for the crossfade.
function getKoppenStripeStyle(zoom) {
const zoomSteps = zoom - DEFAULT_COUNTRY_VIEW.zoom;
// Stripes follow the map: the width doubles with each zoom-in step, so a county's stripe count stays constant.
const followedWidth = getKoppenMixedStripePairWidthPx() * 2 ** zoomSteps;
if (followedWidth < KOPPEN_MIXED_MIN_PAIR_WIDTH_PX) {
// Far-out: double the map-following width (half as many lines) until it reaches the minimum.
let pairWidth = followedWidth * 2;
while (pairWidth < KOPPEN_MIXED_MIN_PAIR_WIDTH_PX) {
pairWidth *= 2;
}
return { name: "far", zoom, pairWidth, primaryFraction: KOPPEN_MIXED_FAR_PRIMARY_STRIPE_FRACTION };
}
return { name: "standard", zoom, pairWidth: followedWidth, primaryFraction: KOPPEN_MIXED_PRIMARY_STRIPE_FRACTION };
}
// Build one diagonal stripe pattern: a primary-color tile with secondary-color bands on top.
function buildKoppenStripePattern(patternId, primaryColor, secondaryColor) {
const pattern = document.createElementNS(SVG_NAMESPACE, "pattern");
pattern.setAttribute("id", patternId);
pattern.setAttribute("patternUnits", "userSpaceOnUse");
pattern.setAttribute("patternTransform", `rotate(${KOPPEN_MIXED_STRIPE_ROTATION_DEGREES})`);
const background = document.createElementNS(SVG_NAMESPACE, "rect");
background.setAttribute("class", "koppen-stripe-primary");
background.setAttribute("x", "0");
background.setAttribute("y", "0");
background.setAttribute("fill", primaryColor);
const secondaries = document.createElementNS(SVG_NAMESPACE, "g");
secondaries.setAttribute("class", "koppen-stripe-secondaries");
secondaries.setAttribute("fill", secondaryColor);
pattern.append(background, secondaries);
return pattern;
}
// Return a style's secondary bands repeated across a tile, as [from, to] pixel ranges.
function getKoppenSecondaryBands(pairWidth, primaryFraction, tileWidth) {
const bands = [];
const pairCount = Math.round(tileWidth / pairWidth);
for (let index = 0; index < pairCount; index += 1) {
bands.push([(index + primaryFraction) * pairWidth, (index + 1) * pairWidth]);
}
return bands;
}
// Split a tile into bands that are secondary in both layouts ("shared"), only the old one, or only the new one.
function classifyKoppenStripeBands(oldBands, newBands, tileWidth) {
const isInside = (bands, x) => bands.some(([from, to]) => x > from && x < to);
const edges = Array.from(new Set([0, tileWidth, ...oldBands.flat(), ...newBands.flat()])).sort((a, b) => a - b);
const segments = [];
for (let index = 0; index < edges.length - 1; index += 1) {
const from = edges[index];
const to = edges[index + 1];
if (to - from < 1e-6) {
continue;
}
const middle = (from + to) / 2;
const isOld = isInside(oldBands, middle);
const isNew = isInside(newBands, middle);
const role = isOld && isNew ? "shared" : isOld ? "old" : isNew ? "new" : null;
if (!role) {
continue;
}
const last = segments[segments.length - 1];
if (last && last.role === role && Math.abs(last.to - from) < 1e-6) {
last.to = to;
} else {
segments.push({ role, from, to });
}
}
return segments;
}
// Draw secondary bands in every stripe pattern, sized to one tile. New-only bands start hidden.
function drawKoppenStripeBands(tileWidth, segments) {
koppenStripeDefs?.querySelectorAll("pattern").forEach((pattern) => {
pattern.setAttribute("width", String(tileWidth));
pattern.setAttribute("height", String(tileWidth));
const background = pattern.querySelector(".koppen-stripe-primary");
background.setAttribute("width", String(tileWidth));
background.setAttribute("height", String(tileWidth));
const rects = segments.map(({ role, from, to }) => {
const rect = document.createElementNS(SVG_NAMESPACE, "rect");
rect.setAttribute("x", String(from));
rect.setAttribute("y", "0");
rect.setAttribute("width", String(to - from));
rect.setAttribute("height", String(tileWidth));
rect.setAttribute("data-role", role);
rect.setAttribute("fill-opacity", role === "new" ? "0" : "1");
return rect;
});
pattern.querySelector(".koppen-stripe-secondaries").replaceChildren(...rects);
});
anchorKoppenStripePatterns();
}
// Anchor every stripe pattern to the map projection's fixed origin rather than Leaflet's pixel
// origin, which moves on every zoom, so the stripes stay fixed to the ground. The offset is reduced
// to within one tile, which draws identical stripes and avoids precision loss at deep zooms.
function anchorKoppenStripePatterns() {
if (!koppenStripeDefs) {
return;
}
// Leaflet's SVG coordinates are projected pixels minus the pixel origin, so the projection's
// origin sits at minus the pixel origin. Express that point in the rotated pattern's coordinates.
const origin = map.getPixelOrigin();
const angle = (-KOPPEN_MIXED_STRIPE_ROTATION_DEGREES * Math.PI) / 180;
const offsetX = -origin.x * Math.cos(angle) + origin.y * Math.sin(angle);
const offsetY = -origin.x * Math.sin(angle) - origin.y * Math.cos(angle);
koppenStripeDefs.querySelectorAll("pattern").forEach((pattern) => {
const tileWidth = Number(pattern.getAttribute("width"));
if (!(tileWidth > 0)) {
return;
}
const shiftX = ((offsetX % tileWidth) + tileWidth) % tileWidth;
const shiftY = ((offsetY % tileWidth) + tileWidth) % tileWidth;
pattern.setAttribute(
"patternTransform",
`rotate(${KOPPEN_MIXED_STRIPE_ROTATION_DEGREES}) translate(${shiftX} ${shiftY})`
);
});
}
// Draw a stripe style's final layout: one secondary band per stripe pair.
function applyKoppenStripeStyle(style) {
const bandStart = style.pairWidth * style.primaryFraction;
drawKoppenStripeBands(style.pairWidth, [{ role: "shared", from: bandStart, to: style.pairWidth }]);
koppenStripeStyle = style;
}
// Crossfade from the old stripes, exactly as the zoom animation left them, to a new style.
// Bands secondary in both layouts stay solid; old-only bands fade out while new-only bands fade in.
function crossfadeKoppenStripes(oldPairWidth, oldPrimaryFraction, style) {
const tileWidth = Math.max(oldPairWidth, style.pairWidth);
const segments = classifyKoppenStripeBands(
getKoppenSecondaryBands(oldPairWidth, oldPrimaryFraction, tileWidth),
getKoppenSecondaryBands(style.pairWidth, style.primaryFraction, tileWidth),
tileWidth
);
drawKoppenStripeBands(tileWidth, segments);
koppenStripeStyle = style;
const oldBands = koppenStripeDefs.querySelectorAll('rect[data-role="old"]');
const newBands = koppenStripeDefs.querySelectorAll('rect[data-role="new"]');
const fadeStart = performance.now();
const fadeStep = (now) => {
// The frame timestamp can fall slightly before fadeStart, so clamp to 0..1.
const progress = Math.min(Math.max((now - fadeStart) / KOPPEN_MIXED_STYLE_FADE_MS, 0), 1);
oldBands.forEach((rect) => rect.setAttribute("fill-opacity", String(1 - progress)));
newBands.forEach((rect) => rect.setAttribute("fill-opacity", String(progress)));
if (progress < 1) {
koppenStripeFadeFrame = requestAnimationFrame(fadeStep);
} else {
koppenStripeFadeFrame = null;
applyKoppenStripeStyle(style);
}
};
koppenStripeFadeFrame = requestAnimationFrame(fadeStep);
}
// True when two widths differ by a whole power of two, so both layouts fit one tile.
function canCrossfadeKoppenStripes(widthA, widthB) {
const exponent = Math.log2(Math.max(widthA, widthB) / Math.min(widthA, widthB));
return Math.abs(exponent - Math.round(exponent)) < 0.01;
}
// Resize stripes after a zoom, crossfading from the old layout when the stripe size changes.
function updateKoppenStripesForZoom() {
if (!koppenStripeDefs) {
return;
}
const zoom = map.getZoom();
const style = getKoppenStripeStyle(zoom);
const previous = koppenStripeStyle;
if (koppenStripeFadeFrame !== null) {
cancelAnimationFrame(koppenStripeFadeFrame);
koppenStripeFadeFrame = null;
}
if (!previous) {
applyKoppenStripeStyle(style);
return;
}
// Leaflet's zoom animation stretched the old stripes by 2 to the power of the zoom change.
const oldPairWidth = previous.pairWidth * 2 ** (zoom - previous.zoom);
// Fade only when the stripe size changes; otherwise the new layout is drawn instantly.
const sizeChanges = Math.abs(oldPairWidth - style.pairWidth) >= 0.01;
if (!sizeChanges || !canCrossfadeKoppenStripes(oldPairWidth, style.pairWidth)) {
applyKoppenStripeStyle(style);
return;
}
crossfadeKoppenStripes(oldPairWidth, previous.primaryFraction, style);
}
// Create one stripe pattern for each primary/secondary pair used by Mixed counties.
function createKoppenStripePatterns() {
if (!koppenStripeDefs) {
// Patterns live in a hidden SVG; url(#id) fills resolve them anywhere in the page.
const svg = document.createElementNS(SVG_NAMESPACE, "svg");
svg.setAttribute("width", "0");
svg.setAttribute("height", "0");
svg.setAttribute("aria-hidden", "true");
svg.style.position = "absolute";
koppenStripeDefs = document.createElementNS(SVG_NAMESPACE, "defs");
svg.appendChild(koppenStripeDefs);
document.body.appendChild(svg);
map.on("zoomend", updateKoppenStripesForZoom);
// Leaflet can also move its pixel origin without a zoom (after long pans or re-centring).
map.on("viewreset moveend", anchorKoppenStripePatterns);
}
const patterns = new Map();
countyStats.forEach((stats) => {
const classes = getKoppenStripeClasses(stats);
if (!classes) {
return;
}
const patternId = getKoppenStripePatternId(classes.primary, classes.secondary);
if (!patterns.has(patternId)) {
const primaryColor = KOPPEN_CLASS_META[classes.primary].color;
const secondaryColor = KOPPEN_CLASS_META[classes.secondary].color;
patterns.set(patternId, buildKoppenStripePattern(patternId, primaryColor, secondaryColor));
}
});
koppenStripeDefs.replaceChildren(...patterns.values());
applyKoppenStripeStyle(getKoppenStripeStyle(map.getZoom()));
}
// Pick a county's map fill: a stripe pattern for Mixed Koppen counties, otherwise the value color.
function fillForCounty(stats) {
if (currentMetricKey === "koppenZone") {
const classes = getKoppenStripeClasses(stats);
if (classes) {
return `url(#${getKoppenStripePatternId(classes.primary, classes.secondary)})`;
}
}
return colorForValue(currentMetricKey, stats[currentMetricKey]);
}
// Build the Mixed Climate legend swatch: gray stripes in the same split as the map.
function buildKoppenMixedLegendSwatchBackground() {
const periodPx = 6;
const primaryPx = periodPx * KOPPEN_MIXED_PRIMARY_STRIPE_FRACTION;
return (
`repeating-linear-gradient(45deg, ${KOPPEN_MIXED_LEGEND_PRIMARY_GRAY} 0 ${primaryPx}px, ` +
`${KOPPEN_MIXED_LEGEND_SECONDARY_GRAY} ${primaryPx}px ${periodPx}px)`
);
}
// Format a county's metric value, naming the stripe classes for Mixed Koppen counties.
function formatCountyMetricValue(metricKey, stats) {
const text = formatMetricValue(metricKey, stats[metricKey]);
const classes = metricKey === "koppenZone" ? getKoppenStripeClasses(stats) : null;
return classes ? `${text} (${classes.primary} primary, ${classes.secondary} secondary)` : text;
}
// Check whether a county should appear active under the current filter. // Check whether a county should appear active under the current filter.
function shouldFeaturePassFilter(stats) { function shouldFeaturePassFilter(stats) {
if (isCurrentMetricCategorical()) { if (isCurrentMetricCategorical()) {
@@ -1113,7 +1456,12 @@ function shouldFeaturePassFilter(stats) {
if (!categoryValue) { if (!categoryValue) {
return false; return false;
} }
return currentCategoricalFilter === "all" || categoryValue === currentCategoricalFilter; if (currentCategoricalFilter === "all" || categoryValue === currentCategoricalFilter) {
return true;
}
// A Mixed Koppen county also matches its primary and secondary classes.
const stripeClasses = currentMetricKey === "koppenZone" ? getKoppenStripeClasses(stats) : null;
return Boolean(stripeClasses) && [stripeClasses.primary, stripeClasses.secondary].includes(currentCategoricalFilter);
} }
const rawValue = stats[currentMetricKey]; const rawValue = stats[currentMetricKey];
@@ -1141,13 +1489,8 @@ function styleForCounty(feature) {
} }
const inRange = shouldFeaturePassFilter(stats); const inRange = shouldFeaturePassFilter(stats);
const metricValue = stats[currentMetricKey];
return buildCountyStyle( return buildCountyStyle(isSelected, inRange ? fillForCounty(stats) : NO_DATA_FILL_COLOR, inRange ? 0.72 : 0.25);
isSelected,
inRange ? colorForValue(currentMetricKey, metricValue) : NO_DATA_FILL_COLOR,
inRange ? 0.72 : 0.25
);
} }
// Apply shared stroke styling around a county fill style. // Apply shared stroke styling around a county fill style.
@@ -1171,10 +1514,13 @@ function updateLegend() {
const rows = values const rows = values
.map((value) => { .map((value) => {
const selectedClass = currentCategoricalFilter === value ? " legend-koppen-row-selected" : ""; const selectedClass = currentCategoricalFilter === value ? " legend-koppen-row-selected" : "";
const isMixed = currentMetricKey === "koppenZone" && value === KOPPEN_MIXED_CLASS;
const swatchBackground = isMixed ? buildKoppenMixedLegendSwatchBackground() : colorForValue(currentMetricKey, value);
const mixedNote = isMixed ? "<br><small>Stripes show each county's top two classes</small>" : "";
return ` return `
<div class="legend-koppen-row${selectedClass}"> <div class="legend-koppen-row${selectedClass}">
<span class="legend-koppen-dot" style="background:${colorForValue(currentMetricKey, value)}"></span> <span class="legend-koppen-dot" style="background:${swatchBackground}"></span>
<span>${getCategoricalDisplayLabel(currentMetricKey, value)}</span> <span>${getCategoricalDisplayLabel(currentMetricKey, value)}${mixedNote}</span>
</div> </div>
`; `;
}) })
@@ -1219,7 +1565,7 @@ function buildCountyDetailsBody(meta, stats) {
// Build one metric row for the selected-county details panel. // Build one metric row for the selected-county details panel.
function buildMetricDetailRow(metricKey, stats) { function buildMetricDetailRow(metricKey, stats) {
return `<p><strong>${getMetricConfig(metricKey).label}:</strong> ${formatMetricValue(metricKey, stats[metricKey])}</p>`; return `<p><strong>${getMetricConfig(metricKey).label}:</strong> ${formatCountyMetricValue(metricKey, stats)}</p>`;
} }
// Update the selected-county panel for the current selection state. // Update the selected-county panel for the current selection state.
@@ -1433,6 +1779,30 @@ function updateRangeLabels() {
maxValue.textContent = formatMetricValue(currentMetricKey, max); maxValue.textContent = formatMetricValue(currentMetricKey, max);
} }
// Move a hovered range input by one configured step without scrolling the page.
function adjustRangeOnWheel(event) {
if (event.deltaY === 0 || isCurrentMetricCategorical()) {
return;
}
const input = event.currentTarget;
const previousValue = input.value;
input.focus({ preventScroll: true });
if (event.deltaY < 0) {
input.stepUp();
} else {
input.stepDown();
}
if (input.value === previousValue) {
return;
}
event.preventDefault();
input.dispatchEvent(new Event("input", { bubbles: true }));
}
// Zoom the map to a county without animation. // Zoom the map to a county without animation.
function focusCounty(layer) { function focusCounty(layer) {
map.stop(); map.stop();
@@ -1454,9 +1824,9 @@ function openCountyPopup(layer, countyFips) {
layer layer
.bindPopup( .bindPopup(
`<strong>${meta.displayName}</strong><br>${METRICS[currentMetricKey].label}: ${formatMetricValue( `<strong>${meta.displayName}</strong><br>${METRICS[currentMetricKey].label}: ${formatCountyMetricValue(
currentMetricKey, currentMetricKey,
stats[currentMetricKey] stats
)}`, )}`,
{ autoPan: false } { autoPan: false }
) )
@@ -1636,7 +2006,7 @@ function configureRangeControls() {
function configureRangeInput(input, metric, value) { function configureRangeInput(input, metric, value) {
input.min = String(metric.min); input.min = String(metric.min);
input.max = String(metric.max); input.max = String(metric.max);
input.step = String(metric.sliderStep ?? 1); input.step = String(getMetricSliderStep(metric));
input.value = String(value); input.value = String(value);
} }
@@ -1789,6 +2159,7 @@ async function loadCounties() {
const geojson = await response.json(); const geojson = await response.json();
geojson.features = geojson.features.filter(prepareCountyFeature); geojson.features = geojson.features.filter(prepareCountyFeature);
createKoppenStripePatterns();
applyNumericMetricBoundsFromData(); applyNumericMetricBoundsFromData();
configureCategoricalFilterControl(); configureCategoricalFilterControl();
@@ -1845,6 +2216,8 @@ function configureDomEventListeners() {
metricSelect.addEventListener("change", onMetricChange); metricSelect.addEventListener("change", onMetricChange);
minRange.addEventListener("input", restyleAllCounties); minRange.addEventListener("input", restyleAllCounties);
maxRange.addEventListener("input", restyleAllCounties); maxRange.addEventListener("input", restyleAllCounties);
minRange.addEventListener("wheel", adjustRangeOnWheel, { passive: false });
maxRange.addEventListener("wheel", adjustRangeOnWheel, { passive: false });
resetViewButton.addEventListener("click", resetMapView); resetViewButton.addEventListener("click", resetMapView);
clearSelectionButton.addEventListener("click", clearSelection); clearSelectionButton.addEventListener("click", clearSelection);
+3222 -3222
View File
File diff suppressed because it is too large Load Diff
+620
View File
@@ -0,0 +1,620 @@
# County Filter Calculations
This document records the equations, implementation behavior, and review notes
for the 12 county filters currently exposed by the climate explorer.
## Scope and notation
The current metric list is defined in `app.js` under `METRICS`. Unless noted
otherwise, long-term climate metrics use the 1991--2020 reference period.
| Symbol | Meaning |
| --- | --- |
| \(c\) | County |
| \(i\) | Raster or model grid cell |
| \(m\) | Calendar month |
| \(d\) | Calendar day |
| \(h\) | Hour or NSRDB time row |
| \(y\) | Year |
| \(G_c\) | Valid grid cells assigned to county \(c\) |
| \(V_c\) | Valid observations for county \(c\) |
| \(\mathbf{1}[A]\) | 1 when condition \(A\) is true; otherwise 0 |
| \(\operatorname{clip}(x,a,b)\) | Restrict \(x\) to the interval \([a,b]\) |
Missing values are omitted from means unless a metric-specific rule below says
otherwise.
## Browser filter predicate
For every numeric metric, a county is active when its value is finite and falls
inside the selected range, including both endpoints:
$$
\operatorname{passes}(c) \iff L \le x_c \le U.
$$
For a categorical metric:
$$
\operatorname{passes}(c) \iff
\left(F=\text{All}\right) \lor \left(x_c=F\right).
$$
The Köppen-Geiger filter also matches Mixed counties by their top two classes;
see §1.
Counties with missing categorical values, null numeric values, or non-finite
numeric values do not pass. The numeric slider limits are derived from the
loaded data and rounded outward by each metric's configured `boundsStep`.
Implementation: `shouldFeaturePassFilter` and `getUniqueCategoryValuesInData`
in `app.js`.
## Filter inventory
| Group | UI label | Data key | Type | Display unit |
| --- | --- | --- | --- | --- |
| Climate Classification | Köppen-Geiger Climate Class | `koppenZone` | Categorical | Class |
| Temperature & Extremes | Annual Avg Temperature (Normals) | `avgTempF` | Numeric | °F |
| Temperature & Extremes | Diurnal Temperature Range | `avgDiurnalTempRangeF` | Numeric | °F difference |
| Temperature & Extremes | Annual Extreme Temperature Days | `absoluteExtremeDays` | Numeric | Days/year |
| Temperature & Extremes | Annual 90 F+ Heat Index Days | `humidHeatDays` | Numeric | Days/year |
| Precipitation & Moisture | Annual Precipitation (Normals) | `annualPrecipIn` | Numeric | Inches/year |
| Precipitation & Moisture | Seasonality Index | `seasonalityIndex` | Numeric | 0--100 index |
| Precipitation & Moisture | Wettest Month | `wettestPrecipMonth` | Categorical | Month |
| Precipitation & Moisture | Driest Month | `driestPrecipMonth` | Categorical | Month |
| Precipitation & Moisture | Summer Specific Humidity | `avgSummerSpecificHumidityGKg` | Numeric | g/kg |
| Solar Resource | Mean Daily Global Horizontal Radiation (GHI) | `meanDailyGlobalHorizontalRadiationKwhM2Day` | Numeric | kWh/m²/day |
| Solar Resource | Clear-Sky GHI Reduction Index | `clearSkyGhiReductionIndex` | Numeric | Ratio |
## 1. Köppen-Geiger Climate Class
**Data keys:** `koppenZone`; `koppenPrimaryClass` and `koppenSecondaryClass`
for Mixed counties
There is no continuous numerical score. Each county is classified from the
share of its land covered by each Köppen class. The rule was adopted on
2026-09-12 and applied to `data/climate-data.csv` on 2026-09-13.
Let \(s_{c,k}\) be the share of county \(c\)'s land area covered by class \(k\).
Each valid raster cell is weighted by the area of the cell that lies inside the
county, \(a_{c,i}\); ocean and no-data cells are excluded:
$$
s_{c,k}
=
\frac{\sum_{i}a_{c,i}\,\mathbf{1}[K_i=k]}{\sum_{i}a_{c,i}}.
$$
Rank the classes so that \(s_{c,(1)}\ge s_{c,(2)}\ge\cdots\), with
\(s_{c,(2)}=0\) when only one class is present. A county is **predominantly**
class \(k_{(1)}\), shown as a single color, if and only if both conditions hold:
$$
K_c=
\begin{cases}
k_{(1)}, & s_{c,(1)}\ge0.50 \;\text{and}\; s_{c,(1)}-s_{c,(2)}\ge0.05,\\
\text{Mixed}, & \text{otherwise}.
\end{cases}
$$
The gap is measured in percentage points. A county that fails either condition
is classified as **Mixed** (shown as "Mixed Climate"). For Mixed counties,
\(p_c=k_{(1)}\) and \(q_c=k_{(2)}\) are stored in `koppenPrimaryClass` and
`koppenSecondaryClass`, and the map draws the county with stripes of those two
classes (see [koppen-mixed-display-plan.md](koppen-mixed-display-plan.md)).
Both columns are blank for predominant counties. A county with no valid raster
cells is left blank; none in the 50 states and DC is.
Each county is read from a small raster window around its polygon. Each cell is
split into 16 × 16 sub-cells to estimate the fraction inside the county, and
scaled by \(\cos(\text{latitude})\) for its true surface area. A county whose
polygon crosses the 180th meridian (Aleutians West, AK) is split into one piece
on each side, and each piece is read from its own window.
**Filter.** Choosing a class \(F\) shows counties that are predominantly that
class and Mixed counties where it is the primary or secondary class; the
"Mixed Climate" option shows every Mixed county:
$$
\operatorname{passes}(c) \iff
\left(F=\text{All}\right) \lor \left(K_c=F\right) \lor
\left(K_c=\text{Mixed} \land F\in\{p_c,q_c\}\right).
$$
Implementation: `scripts/build_county_koppen_metric.py` (`rank_class_shares`,
`classify`, `build_koppen_records`) writes `data/metrics/koppen.csv`; area
weighting, raster windows, and the 180th-meridian split are in
`scripts/common/county_zonal_stats.py`;
`scripts/apply_koppen_metric_to_climate_data.py` copies the three columns into
`data/climate-data.csv`; the filter is `shouldFeaturePassFilter` in `app.js`.
**Rationale.** No published standard defines when an area is predominantly one
Köppen class; the classification is defined per grid cell. The 50% condition
means the label describes more than half of the county's land. The 5-point gap
condition catches near 50/50 splits, such as Schenectady, NY (Dfb 50.1%,
Dfa 49.9%), where a single label would rest on a margin of a few tenths of a
point. Map-unit purity standards from other fields were considered and
rejected: the FAO Land Cover Classification System treats a unit as single
only above 80%, and USDA soil survey consociations allow roughly 15--25%
dissimilar inclusions. Applied to counties, those thresholds would mark about
40--47% of the map area as mixed. A published county-level Köppen dataset
(Audirac, Harvard Dataverse, 2024) uses the plurality class and reports the
share of every class, without a threshold.
**Results.** Using the Beck et al. 2023 1991--2020 1 km raster, for the 3,143
counties in the 50 states and DC:
| Classification | Counties | Share of counties | Share of map area |
| --- | --- | --- | --- |
| Predominant (single class) | 3,010 | 95.8% | 85.6% |
| Mixed climate | 133 | 4.2% | 14.4% |
Of the 133 mixed counties, 111 have no class covering 50% or more (37 of
these also have a gap under 5 points), and 22 have a majority class whose
runner-up is within 5 points. They are concentrated in the mountain West and
Alaska: Colorado and Alaska (13 each), California and Montana (12 each),
Washington (10), Idaho (9), and Utah (8). No county's top three classes fall
within 2 points of one another, so no additional mixed categories are needed.
**Compared with the previous method** (below). Every county whose label changed
became Mixed; no county moved to a different single class. The touched-cell
counts and the area-weighted shares pick a different plurality winner in only
3 counties (Denver, CO; Schenectady, NY; Hood River, OR), all of which are
Mixed under this rule. Applying this rule to touched-cell counts instead of
area-weighted shares would classify 9 counties differently, because
touched-cell counts give full weight to boundary cells that are mostly
outside the county.
**Boundary cases.** Five counties lie within 0.25 points of a cutoff that
decides their outcome: Sanpete, UT (Dfb 49.93%, Mixed), Giles, VA (49.95%,
Mixed), Placer, CA (Csa 50.22%), Albany, WY (Dfb 50.24%), and Park, MT (gap
5.23 points, Dfb). Their classification depends on the precision of the
area weighting. Aleutians West, AK (02016) spans the antimeridian, so its
shares were computed without sub-cell sampling. Puerto Rico is not covered by
these figures; it has 9 more Mixed counties and is not shown in the app.
### Previous method
Until 2026-09-13, `koppenZone` was the most frequent valid Köppen raster code
among the cells touched by the county polygon:
$$
K_c = \underset{k}{\arg\max}\;
\sum_{i\in G_c}\mathbf{1}[K_i=k].
$$
Cells were equally weighted, regardless of how much of each cell lay inside the
county; an exact tie went to the smallest raster code, and a county with no
valid cell was assigned `Cfa`. `scripts/build_county_climate_data.py` still
computes this value until that script is retired (pipeline plan, Phase 3), so
the Köppen apply step must run after it.
## 2. Annual Avg Temperature (Normals)
**Data key:** `avgTempF`
If the NOAA source is a historical monthly series, the script first forms a
1991--2020 climatology for each calendar month and cell:
$$
T_{i,m}^{\mathrm{norm}}
=
\frac{1}{Y_{i,m}}
\sum_{y\in V_{i,m}}T_{i,m,y}.
$$
The monthly county value is an unweighted mean of valid touched raster cells:
$$
T_{c,m}
=
\frac{1}{|G_{c,m}|}
\sum_{i\in G_{c,m}}T_{i,m}^{\mathrm{norm}}.
$$
The annual value is the equally weighted mean of the available monthly values,
converted from Celsius to Fahrenheit:
$$
T_c(^\circ\mathrm{F})
=
\left(
\frac{1}{M_c}\sum_{m\in V_c}T_{c,m}
\right)\frac{9}{5}+32.
$$
Normally \(M_c=12\). The stored value is rounded to 0.1 °F.
Implementation: `scripts/build_county_climate_data.py:124--166`,
`scripts/build_county_climate_data.py:206--221`, and
`scripts/build_county_climate_data.py:478--514`.
## 3. Diurnal Temperature Range
**Data key:** `avgDiurnalTempRangeF`
For each day with paired county Tmax and Tmin:
$$
DTR_{c,d}=T^{\max}_{c,d}-T^{\min}_{c,d}.
$$
Days with a missing input or a negative range are excluded. The final metric is
the mean across all retained days in 1991--2020, followed by conversion of a
Celsius temperature *difference* to a Fahrenheit difference:
$$
\overline{DTR}_c(^\circ\mathrm{F})
=
\frac{9}{5}
\left(
\frac{1}{N_c}\sum_{d\in V_c}DTR_{c,d}
\right).
$$
There is correctly no \(+32\) term when converting a temperature difference.
The output artifact stores two decimal places.
Implementation: `scripts/build_county_diurnal_temperature_range.py:46--104`
and `scripts/build_county_diurnal_temperature_range.py:132--151`.
## 4. Annual Extreme Temperature Days
**Data key:** `absoluteExtremeDays`
A valid paired day is counted once when either the hot or cold absolute threshold
is met:
$$
E_{c,y}
=
\sum_{d\in V_{c,y}}
\mathbf{1}\!\left[
T^{\max}_{c,d}\ge95^\circ\mathrm{F}
\;\lor\;
T^{\min}_{c,d}\le0^\circ\mathrm{F}
\right].
$$
The county filter is the arithmetic mean of the yearly counts:
$$
E_c=\frac{1}{Y_c}\sum_{y\in V_c}E_{c,y}.
$$
The checked-in data uses 1991--2025. The Boolean OR means that a hypothetical
day meeting both conditions is still counted only once. A year is included when
an annual record exists; counts are not normalized to 365 or 366 valid days.
Implementation: `scripts/build_county_locally_extreme_data.py:473--540`,
`scripts/build_county_locally_extreme_data.py:669--674`, and
`scripts/build_county_locally_extreme_data.py:705--722`.
## 5. Annual 90 F+ Heat Index Days
**Data key:** `humidHeatDays`
The daily proxy pairs NOAA nClimGrid-Daily county Tmax, \(T\), with the gridMET
county daily minimum relative humidity, \(R\). Relative humidity is clipped to
\([0,100]\).
The NWS simple Heat Index estimate is calculated in two steps:
$$
S=\frac{1}{2}\left[T+61+1.2(T-68)+0.094R\right],
$$
$$
HI_s=\frac{S+T}{2}.
$$
When \(HI_s\ge80^\circ\mathrm{F}\), the Rothfusz regression is used:
$$
\begin{aligned}
HI_r={}&-42.379+2.04901523T+10.14333127R-0.22475541TR\\
&-0.00683783T^2-0.05481717R^2+0.00122874T^2R\\
&+0.00085282TR^2-0.00000199T^2R^2.
\end{aligned}
$$
For \(R<13\) and \(80\le T\le112\), subtract:
$$
A_{low}
=
\frac{13-R}{4}
\sqrt{\max\left(\frac{17-|T-95|}{17},0\right)}.
$$
For \(R>85\) and \(80\le T\le87\), add:
$$
A_{high}=\frac{R-85}{10}\frac{87-T}{5}.
$$
Thus:
$$
HI(T,R)=
\begin{cases}
HI_s, & HI_s<80,\\
HI_r-A_{low}, & HI_s\ge80 \text{ and the low-RH condition holds},\\
HI_r+A_{high}, & HI_s\ge80 \text{ and the high-RH condition holds},\\
HI_r, & \text{otherwise}.
\end{cases}
$$
The yearly and long-term values are:
$$
H_{c,y}=\sum_{d\in V_{c,y}}
\mathbf{1}[HI(T^{\max}_{c,d},R^{\min}_{c,d})\ge90],
$$
$$
H_c=\frac{1}{Y_c}\sum_{y\in V_c}H_{c,y}.
$$
The period is 1991--2020. A year with at least one valid paired day contributes
equally to the final average; there is no completeness adjustment.
Implementation: `scripts/summarize_county_gridmet_humidity.py:439--474` and
`scripts/summarize_county_gridmet_humidity.py:547--617`.
## 6. Annual Precipitation (Normals)
**Data key:** `annualPrecipIn`
Each monthly county total, \(P_{c,m}\), is the unweighted mean of valid raster
cells touched by the county. Annual precipitation is the sum of available
monthly totals, converted from millimeters to inches:
$$
P_c(\mathrm{in})
=
\frac{1}{25.4}
\sum_{m\in V_c}P_{c,m}(\mathrm{mm}).
$$
The stored value is rounded to 0.1 inch. The implementation only requires one
valid month, so missing months produce a partial annual sum rather than a blank.
Implementation: `scripts/build_county_climate_data.py:401--412` and
`scripts/build_county_climate_data.py:478--511`.
## 7. Seasonality Index
**Data key:** `seasonalityIndex`
The filter is the coefficient of variation of available monthly precipitation
totals. First calculate the monthly mean and population standard deviation:
$$
\mu_c=\frac{1}{M_c}\sum_{m\in V_c}P_{c,m},
$$
$$
\sigma_c=
\sqrt{
\frac{1}{M_c}
\sum_{m\in V_c}(P_{c,m}-\mu_c)^2
}.
$$
Then:
$$
SI_c=
\operatorname{round}\!\left(
\operatorname{clip}\!\left(
100\frac{\sigma_c}{\mu_c},0,100
\right)
\right).
$$
If \(\mu_c\le0\), the index is set to zero. The value is stored as an integer.
Implementation: `scripts/build_county_climate_data.py:487--521`.
## 8. Wettest Month
**Data key:** `wettestPrecipMonth`
$$
W_c=\underset{m\in V_c}{\arg\max}\;P_{c,m}.
$$
Missing monthly values are ignored. An exact tie resolves to the earliest tied
month because `numpy.nanargmax` returns the first occurrence.
Implementation:
`scripts/apply_precipitation_month_metrics_to_climate_data.py:56--90`.
## 9. Driest Month
**Data key:** `driestPrecipMonth`
$$
D_c=\underset{m\in V_c}{\arg\min}\;P_{c,m}.
$$
Missing monthly values are ignored. An exact tie likewise resolves to the
earliest tied month.
Implementation:
`scripts/apply_precipitation_month_metrics_to_climate_data.py:56--90`.
## 10. Summer Specific Humidity
**Data key:** `avgSummerSpecificHumidityGKg`
gridMET cells whose centers fall inside a county are weighted by the cosine of
their latitude to approximate their relative surface areas on a latitude--longitude
grid:
$$
q_{c,d}
=
\frac{
\sum_{i\in G_{c,d}}q_{i,d}\cos(\phi_i)
}{
\sum_{i\in G_{c,d}}\cos(\phi_i)
}.
$$
The metric averages the valid daily county values for June, July, and August,
then converts kg/kg to g/kg:
$$
q_c^{\mathrm{summer}}
=
1000\left(
\frac{1}{N_c}
\sum_{d\in V_c,\;m(d)\in\{6,7,8\}}q_{c,d}
\right).
$$
If no grid-cell center falls inside a county, the nearest grid cell to an
interior representative point is used.
Implementation: `scripts/summarize_county_gridmet_humidity.py:297--380`,
`scripts/summarize_county_gridmet_humidity.py:409--434`, and
`scripts/summarize_county_gridmet_humidity.py:555--605`.
## 11. Mean Daily Global Horizontal Radiation (GHI)
**Data key:** `meanDailyGlobalHorizontalRadiationKwhM2Day`
For each NSRDB site, the current 60-minute, non-leap-year TMY data is summarized
as:
$$
G_s
=
\frac{\sum_h GHI_{s,h}}{1000\times365}
\quad\mathrm{kWh/m^2/day}.
$$
For a county with polygon archive coverage:
$$
G_c
=
\frac{\sum_s A_{c,s}G_s}{\sum_s A_{c,s}},
$$
where \(A_{c,s}\) is the estimated overlap area between the county geometry and
the 4 km square grid cell centered on site \(s\). A representative-point value
is used when a polygon summary is unavailable.
Implementation: `scripts/summarize_nsrdb_county_polygon_archives.py:189--192`
and `scripts/summarize_nsrdb_county_polygon_archives.py:335--378`.
## 12. Clear-Sky GHI Reduction Index
**Data key:** `clearSkyGhiReductionIndex`
Only rows with valid observed and clear-sky GHI and
\(CSGHI_{s,h}\ge50\;\mathrm{W/m^2}\) are treated as daylight rows. For each
retained row:
$$
r_{s,h}
=
\operatorname{clip}\!\left(
\frac{GHI_{s,h}}{CSGHI_{s,h}},0,1
\right).
$$
The site reduction index is:
$$
R_s=1-\frac{1}{N_s}\sum_{h\in V_s}r_{s,h}.
$$
The polygon county value is overlap-area-weighted:
$$
R_c=\frac{\sum_s A_{c,s}R_s}{\sum_s A_{c,s}}.
$$
A representative-point index is used where a polygon summary is unavailable.
This definition is the mean of time-row ratios; it is not generally equal to
\(1-\sum GHI/\sum CSGHI\).
Implementation:
`scripts/summarize_nsrdb_county_polygon_cloud_archives.py:138--218` and
`scripts/summarize_nsrdb_county_polygon_cloud_archives.py:222--304`.
## Calculation review findings
1. **Annual temperature weights months equally.** February has the same weight
as January or July. If the intended label means an average across all days,
monthly normals should instead be weighted by the number of days in each
month.
2. **Base NOAA aggregation is not area-weighted.** Every touched raster cell
receives equal weight, including cells that intersect only a small portion
of a county. This can matter most for small or narrow counties and along
coastlines. *Resolved for Köppen on 2026-09-13: class shares are now
area-weighted (§1).*
3. **Partial precipitation years are accepted.** One valid monthly precipitation
value is sufficient to produce `annualPrecipIn`; absent months silently lower
the annual sum. Requiring all 12 months, or recording completeness, would be
safer.
4. **The Köppen fallback can create false data.** A county with no valid raster
cells is labeled `Cfa` instead of missing. A null value plus an audit flag
would distinguish missing coverage from a genuine humid-subtropical class.
*Resolved on 2026-09-13: the Köppen builder leaves such counties blank. The
fallback remains only in `build_county_climate_data.py`, whose Köppen
value is replaced by the apply step.*
5. **Extreme-day counts are not completeness-normalized.** A partially observed
year contributes a raw count and receives the same weight as a complete year.
Consider requiring a minimum number of valid days or annualizing partial
counts explicitly.
6. **The absolute-extreme metric depends on unrelated percentile thresholds.**
`build_annual_counts` skips a county when its retired local p95/p05 thresholds
are missing, even though the active 95 °F / 0 °F calculation does not require
those percentiles. The absolute calculation should be separated from that
prerequisite.
7. **Heat Index days are a daily-extrema proxy.** Daily Tmax and daily minimum
relative humidity are paired even though their observation times may differ.
The result should not be described as an observed hourly maximum Heat Index.
8. **Heat-year completeness is permissive.** Any year with at least one valid
Tmax/RH pair is included in the equal-year average. A minimum valid-day rule
would reduce low-biased partial-year counts.
9. **The GHI formula assumes hourly, 365-day input.** It is correct for the
current 60-minute, `leap_day=false` requests. If the request interval changes,
the energy sum needs an interval-hours multiplier; leap-day handling would
also need to change the divisor.
10. **Clear-sky reduction averages ratios rather than energy totals.** This is a
valid but specific definition. It gives each retained time row equal weight,
rather than weighting rows by available clear-sky energy. The label and
documentation should retain this distinction.
11. **Spatial weighting is inconsistent across metric families.** Base NOAA
normals use equal touched-cell weights, Köppen uses area-weighted class
shares, gridMET humidity uses
\(\cos(\phi)\) weights on cell centers, and NSRDB polygon metrics use
estimated overlap areas. Cross-metric comparisons should account for these
different county aggregation methods.
## Verification status
This reference was derived from the checked-in calculation and merge scripts,
not solely from UI descriptions. No calculation code was changed. The automated
test suite was not executed during this review because `pytest` is not installed
in either the system Python environment or the project virtual environment.
Section 1 and findings 2, 4, and 11 were updated on 2026-09-14, after the
Köppen classification was reworked and applied.
+198
View File
@@ -0,0 +1,198 @@
# Mixed Climate Display
**Status (2026-09-14):** Complete. The app draws Mixed counties as stripes, the
pipeline writes the two stripe-class columns, and `data/climate-data.csv` has
been updated (142 counties Mixed, 133 of them in the 50 states and DC). The
display was reviewed in the browser and refined (see the history in section 7),
and the project documentation was updated (section 6).
**Scope:** Show Köppen-Geiger "Mixed" counties as striped on the map, filter them
by their top two classes, and carry the two stripe classes from the pipeline into
the app. The classification rule itself is described in
[filter-calculations.md](filter-calculations.md) §1; overall progress is tracked in
[pipeline-plan.md](pipeline-plan.md).
## 1. Design
- **Only Mixed counties are striped.** A county is Mixed when its top class covers
less than 50% of its land or leads the runner-up by less than 5 percentage
points. Every other county keeps its solid class color.
- **Stripe colors** are the county's top class (primary) and runner-up class
(secondary), using the existing Köppen class colors.
- **Stripes are diagonal at 45°**, from upper left to lower right, and run
unbroken to the county boundary. Neighboring Mixed counties share the same
angle and spacing, so stripes line up across borders. The county border is
drawn on top.
- **Stripes follow the map.** The density is exact at the default view
(`DEFAULT_COUNTRY_VIEW`, the project's single reference scale). From there the
stripe pair width doubles with each zoom-in step, all the way to the map's
deepest zoom (19), so the number of stripes in a county stays the same at
every zoom.
- **Far-out style.** When stripes would drop below 5 pixels per pair (zoom 4 and
below), the map-following width is doubled until it reaches 5 pixels (at
least once, so there are half as many lines), giving 8.53-pixel pairs split
70/30 so the secondary stays visible.
- **Stripes are fixed to the ground.** Each pattern is anchored to the map
projection's fixed origin, so a stripe stays on the same ground at every zoom
from 5 to 19. At the far-out zooms every other stripe stays put while the
crossfade blends the rest.
- **Fades happen only when the stripe size changes** (the switch between the
standard and far-out styles, and zoom steps within the far-out range). The fade
is a 250 ms true crossfade that runs just after Leaflet's zoom animation: the
old stripes, exactly as the animation left them, blend into the new ones with
no jump and no gap.
- **Filtering.** Choosing a class shows counties that are predominantly that
class plus Mixed counties where it is the primary or secondary class. Third and
lower classes do not match. One **"Mixed Climate"** option shows all Mixed
counties; there are no per-pair Mixed options.
- **Every class on the map is listed.** Classes that appear only as stripe
colors (today Dsc and Dwc, both in Alaska) are listed in the dropdown and
legend as regular classes.
- **Legend.** The "Mixed Climate" row has a swatch of two neutral grays in the
standard split (lighter gray primary, darker gray secondary), generated from
the split setting so it always matches the map, and a note that stripes show
each county's top two classes.
- **Deferred:** a donut chart of every class share (like Europa Universalis 5).
### Settings
All stripe settings sit together near the top of `app.js`:
| Setting | Value | Controls |
| --- | --- | --- |
| `KOPPEN_MIXED_PRIMARY_STRIPE_FRACTION` | 0.8 | Primary share of each stripe pair, standard style (zoom 5 and closer) |
| `KOPPEN_MIXED_LINE_PAIRS_PER_100_MILES` | 5 | Stripe pairs per 100 miles, measured across the stripes at the default view |
| `KOPPEN_MIXED_MIN_PAIR_WIDTH_PX` | 5 | Below this pair width, stripes switch to the far-out style |
| `KOPPEN_MIXED_FAR_PRIMARY_STRIPE_FRACTION` | 0.7 | Primary share of each stripe pair, far-out style |
| `KOPPEN_MIXED_STYLE_FADE_MS` | 250 | Length of the crossfade |
| `KOPPEN_MIXED_STRIPE_ROTATION_DEGREES` | -45 | Stripe angle (SVG rotation; -45 runs upper left to lower right) |
| `KOPPEN_MIXED_LEGEND_PRIMARY_GRAY`, `KOPPEN_MIXED_LEGEND_SECONDARY_GRAY` | `#d4d8df`, `#6b7383` | Legend swatch colors |
## 2. Data
`data/climate-data.csv` has two columns directly after `koppenZone`:
| Column | Filled when | Value |
| --- | --- | --- |
| `koppenPrimaryClass` | `koppenZone` is `Mixed` | Top class, e.g. `Csb` |
| `koppenSecondaryClass` | `koppenZone` is `Mixed` | Runner-up class, e.g. `Dsb` |
Both are blank for predominant counties. The values come from `koppenTopClass`
and `koppenSecondClass` in `data/metrics/koppen.csv`.
## 3. How the stripes are drawn
The map uses Leaflet 1.9.4 with its default SVG renderer. Each ordered
primary/secondary pair gets one SVG `<pattern>`, kept in a small hidden SVG added
to the page; fill references such as `url(#id)` resolve anywhere in the page, so
Leaflet's own SVG elements are not modified. The 133 Mixed counties in the
50 states and DC use 50 pairs. Order matters: Csb/Dsb and Dsb/Csb give the
primary share to different colors.
- **Pattern contents.** A primary-color background and a group of
secondary-color bands. Normally the group has one band per stripe pair; during
a crossfade it holds shared, old-only, and new-only bands.
- **Coordinates.** `patternUnits="userSpaceOnUse"`, so every county is filled in
the map's coordinate space and stripes align across borders.
- **Anchoring.** `patternTransform="rotate(-45) translate(x y)"`. Leaflet's SVG
coordinates start from a pixel origin that moves on every zoom, so the
translate places the pattern's origin at the map projection's fixed origin,
reduced to within one tile to avoid precision loss at deep zooms.
- **Fill.** A Mixed county's style sets
`fillColor: "url(#koppen-mixed-<primary>-<secondary>)"`.
**Pair width at the default view.** Leaflet uses Web Mercator, where
$$
\text{meters per pixel} = \frac{156{,}543.03 \times \cos(\text{latitude})}{2^{\text{zoom}}}.
$$
At the default view (latitude 39.5°, zoom 5) that is about 3,775 m per pixel,
so 100 miles (160,934 m) is about 42.6 pixels. At 5 pairs per 100 miles, one
stripe pair is about 8.5 pixels: about 6.8 pixels of primary color and
1.7 pixels of secondary color. The app computes this from the settings.
| Zoom | Style | Pair width | Split |
| --- | --- | --- | --- |
| 2–4 | Far-out | 8.53 px | 70/30 |
| 5 (default) | Standard | 8.53 px | 80/20 |
| 6 | Standard | 17.05 px | 80/20 |
| 7 | Standard | 34.11 px | 80/20 |
| 8 | Standard | 68.2 px | 80/20 |
| 9–19 | Standard | Doubling each step, to 139,704 px at zoom 19 | 80/20 |
**Zooming.** During Leaflet's zoom animation the whole SVG layer is scaled with
the map, so the stripes stretch with it. When the zoom ends, the stripes are
recomputed. Because every width is the reference width times a power of two, the
old stripes (previous width × 2 to the power of the zoom change) and the new
stripes fit one pattern tile. If the size is unchanged, the new layout is drawn
directly; if it changed, bands secondary in both layouts stay solid, old-only
bands fade out, and new-only bands fade in.
## 4. App code (`app.js`)
| Area | Functions and settings |
| --- | --- |
| Stripe settings | The `KOPPEN_MIXED_*` constants near the top of the file (section 1) |
| Mixed class and parsing | `KOPPEN_CLASS_META` (`Mixed`, labelled "Mixed Climate", sorted last), `normalizeKoppenCode`, `sanitizeOverrideRecord` (reads the two stripe-class columns) |
| Stripe styles | `getKoppenMixedStripePairWidthPx`, `getKoppenStripeStyle`, `getKoppenStripeClasses`, `getKoppenStripePatternId` |
| Patterns | `createKoppenStripePatterns`, `buildKoppenStripePattern`, `drawKoppenStripeBands`, `applyKoppenStripeStyle`, `anchorKoppenStripePatterns` |
| Zooming and crossfade | `updateKoppenStripesForZoom`, `crossfadeKoppenStripes`, `getKoppenSecondaryBands`, `classifyKoppenStripeBands`, `canCrossfadeKoppenStripes` |
| Fill, filter, and labels | `fillForCounty`, `shouldFeaturePassFilter`, `getUniqueCategoryValuesInData`, `getCategoricalDisplayLabel`, `formatCountyMetricValue` |
| Legend | `updateLegend`, `buildKoppenMixedLegendSwatchBackground` |
The hover tooltip shows only the county name, and `styles.css` is unchanged.
After changing `app.js`, bump `APP_ASSET_VERSION` and the `app.js?v=` query in
`index.html` so browsers load the new files.
## 5. Pipeline
- **`scripts/build_county_koppen_metric.py`** writes `data/metrics/koppen.csv`
with each county's class, top and runner-up class, and their shares.
- **`scripts/apply_koppen_metric_to_climate_data.py`** writes `koppenZone`,
`koppenPrimaryClass`, and `koppenSecondaryClass` into `data/climate-data.csv`,
adding the two stripe columns after `koppenZone` if missing and filling them
only for Mixed counties. `--dry-run` reports changes without writing.
- **`scripts/check_climate_data.py`** requires the two stripe columns: filled if
and only if `koppenZone` is `Mixed`, with valid and different Köppen codes.
- **Tests:** `tests/test_koppen_metric.py` and `tests/test_check_climate_data.py`.
## 6. Documentation (stage 5, done 2026-09-14)
- `filter-calculations.md` §1 describes the applied rule, the stripe columns,
and the filter; the old method is kept as a short "Previous method" note.
Review findings 2 and 4 are marked resolved for Köppen.
- `README.md` and `scripts/county_data_sources.md` describe the new
classification and list the Köppen build, apply, and check commands.
## 7. History
Decisions that were reversed or refined during the browser review:
- **2026-09-13 — Stripes follow the map.** Stripes were first a fixed width on
screen; the owner wanted the number of stripes in a county not to change with
zoom.
- **2026-09-13 — True crossfade.** The first fade faded the stripes out and back
in, which briefly left no stripes; it was replaced by a crossfade from the old
stripes to the new ones.
- **2026-09-13 — Crossfade after the zoom animation.** Starting zoom-in fades
together with the animation was tried and reverted: the perceived slowness had
been the map's own zoom animation.
- **2026-09-13/14 — Near style removed.** A thicker 70/30 split from zoom 9 was
added, then made instant (fades happen only when the stripe size changes), set
to 80/20 by the owner, and finally removed.
- **2026-09-14 — Stripes fixed to the ground.** Stripes slid across the land on
each zoom because patterns were anchored to Leaflet's per-zoom pixel origin.
- **2026-09-14 — No deepest-zoom limit.** A 60-pixel cap, later a 1,100-pixel
safety limit with halving, was removed: any cap forces stripes to thin or
subdivide on screen, and tests showed no performance cost without one.
- **2026-09-14 — Zoom response setting removed.** The fixed-to-the-ground
anchoring and the crossfade both require stripes to follow the map exactly.
## 8. Puerto Rico
The app already leaves Puerto Rico off the map: `prepareCountyFeature` in
`app.js` drops counties with state FIPS 72, and the dropdown and legend are built
from the counties on the map. The rows remain in `data/climate-data.csv` and
`data/metrics/koppen.csv`. Whether to also remove them from the data files is a
separate decision.
+291
View File
@@ -0,0 +1,291 @@
# Data Pipeline Improvement Plan
This plan describes how the county data pipeline will move from scripts that
edit one shared CSV in place to per-metric outputs assembled into the app CSV.
It is a working reference for the filter-by-filter review. Calculation details
for each filter live in [filter-calculations.md](filter-calculations.md).
**Started:** 2026-09-12
**Guiding decision:** review and fix each of the 12 filters one at a time,
confirm each works on its own, and restructure `data/climate-data.csv` only
after all filters are clean. No large rewrite happens up front.
## 1. Current pipeline
### Where each filter comes from
| Filter | Written into `climate-data.csv` by | Upstream scripts |
| --- | --- | --- |
| Köppen-Geiger class (plus the two stripe-class columns) | `apply_koppen_metric_to_climate_data.py` | `build_county_koppen_metric.py` → `data/metrics/koppen.csv` |
| Annual avg temperature | `build_county_climate_data.py` | — |
| Annual precipitation | `build_county_climate_data.py` | — |
| Seasonality index | `build_county_climate_data.py` | — |
| Wettest / driest month | Base build, then overwritten by `apply_precipitation_month_metrics_to_climate_data.py` | — |
| Diurnal temperature range | `apply_diurnal_temperature_range_to_climate_data.py` | `build_county_diurnal_temperature_range.py` |
| Extreme temperature days | `apply_locally_extreme_metric_to_climate_data.py` | `build_county_locally_extreme_data.py` |
| Summer specific humidity | `apply_gridmet_humidity_metric_to_climate_data.py` | `download_gridmet_data.py` → `summarize_county_gridmet_humidity.py` |
| 90 °F+ heat-index days (plus 2 source-FIPS columns) | `apply_gridmet_humidity_metric_to_climate_data.py` | Same as summer humidity |
| Solar GHI | Base build (optional), then replaced by `apply_locally_extreme_metric_to_climate_data.py` | Point: `build_county_representative_points.py` → `fetch_nsrdb_representative_point_ghi.py`. Polygon: `request_nsrdb_county_polygon_ghi_archives.py` → `download_nsrdb_county_polygon_ghi_archives.py` → `summarize_nsrdb_county_polygon_archives.py` |
| Clear-sky GHI reduction | `apply_nsrdb_cloud_metric_to_climate_data.py` | Point: `fetch_nsrdb_representative_point_cloud_metrics.py`. Polygon: `request_nsrdb_county_polygon_cloud_archives.py` → `download_nsrdb_county_polygon_cloud_archives.py` → `summarize_nsrdb_county_polygon_cloud_archives.py` |
Supporting scripts: `request_nsrdb_county_polygon_archives.py` and
`download_nsrdb_county_polygon_archives.py` are the shared engines behind the
GHI and cloud wrappers; `rebuild_nsrdb_representative_point_ghi_summary.py`
rebuilds the point GHI summary from cache; `check_climate_data.py` validates the
final CSV. Shared helpers live in `scripts/common/` (Phase 2). The base build
still writes an old largest-share `koppenZone`, so the Köppen apply step must
run after it.
### Current full-rebuild order
1. `build_county_climate_data.py`
2. `build_county_koppen_metric.py` → `apply_koppen_metric_to_climate_data.py`
3. `apply_precipitation_month_metrics_to_climate_data.py`
4. `build_county_locally_extreme_data.py` → `apply_locally_extreme_metric_to_climate_data.py`
5. `build_county_diurnal_temperature_range.py` → `apply_diurnal_temperature_range_to_climate_data.py`
6. `summarize_county_gridmet_humidity.py` → `apply_gridmet_humidity_metric_to_climate_data.py`
7. `apply_nsrdb_cloud_metric_to_climate_data.py`
The README's enrichment list starts at step 2 and omits steps 1 and 3.
### Problems
1. **Rerunning a step can destroy data.** The base build writes 12 columns,
including the retired `extremeDays`. Later scripts delete, overwrite, or add
columns until the live CSV has 20. Rerunning the base build drops
9 live columns (`koppenPrimaryClass`, `koppenSecondaryClass`,
`avgDiurnalTempRangeF`, `absoluteExtremeDays`,
`clearSkyGhiReductionIndex`, `avgSummerSpecificHumidityGKg`,
`humidHeatDays`, `humidHeatSourceFips`, `humidHeatFipsAdjustment`),
restores `extremeDays`, and rewrites `koppenZone` with the old
largest-share method.
2. **Order is implicit.** The sequence lives in the README, in
`scripts/county_data_sources.md`, and in each script's assumptions.
3. **Column ownership is unclear.** GHI is finalized by the extreme-temperature
apply script; wettest/driest month are computed in two places.
4. **County aggregation is inconsistent.** NOAA uses touched raster cells,
Köppen uses area-weighted shares (since 2026-09-13), gridMET uses cell
centers with cos(latitude) weights, and NSRDB uses overlap areas (finding 11
in `filter-calculations.md`).
5. **No single entry point or final check.** A new user must piece together
about 20 scripts, several large downloads, and an NSRDB API key.
## 2. Target design
1. **One metric, one file.** Each metric pipeline writes a county-level file
under `data/metrics/`, for example `data/metrics/koppen.csv`, containing
`countyFips`, the app value, and any audit columns for that metric.
2. **One assemble step.** A single script joins the metric files into
`data/climate-data.csv`, using `data/metric_sources.json` for the column list
and per-metric source notes, then runs `check_climate_data.py`.
- Run order no longer matters; rerunning one metric cannot damage others.
- Every column has exactly one owner.
- The per-row `source` column moves into `metric_sources.json`.
- Audit columns stay in the metric files rather than the app CSV.
3. **One shared county-aggregation module.** Area-weighted zonal statistics,
including the 180th-meridian split, used by every raster-based metric.
4. **One runner.** For example
`python scripts/pipeline.py --only koppen --skip-download`, with stages for
fetch, build metrics, assemble, and check. Cached downloads are reused by
default.
## 3. Reproduction tiers
The "Reproducing the data" guide (Phase 4) will be organized by how deep a
user needs to go:
| Tier | What the user does | Needs |
| --- | --- | --- |
| 1. Run the app | `.\serve.ps1` with the committed CSV | Nothing else |
| 2. Reassemble | Rebuild `climate-data.csv` from committed metric files | Python environment only |
| 3. Regenerate one metric | Download one source, rebuild one metric file, reassemble | That metric's source data |
| 4. Full rebuild | Everything | All sources; NSRDB API key; large downloads (the NOAA monthly temperature file alone is about 5 GB) and hours of paced NSRDB requests |
The guide will list each dataset's size, download location, API-key needs, and
approximate run time, and the Python requirements will be pinned.
## 4. Roadmap
### Phase 0 — Groundwork (done)
- [x] `scripts/check_climate_data.py` validates the app CSV (8 checks) with
tests in `tests/test_check_climate_data.py`.
- [x] `data/metric_sources.json` created as an empty skeleton.
- [x] Köppen raster reads use a padded window per county, and polygons that
cross the 180th meridian are split (`split_at_antimeridian` in
`scripts/common/county_zonal_stats.py`); output verified identical for
all 3,221 counties; tests in `tests/test_koppen_antimeridian.py`.
### Phase 1 — Filter-by-filter review (in progress)
Each filter goes through the checklist in section 5. Each fix delivers that
metric's own file in `data/metrics/` plus a single-column apply step, so the
existing CSV keeps working until Phase 3.
### Phase 2 — Shared helpers in `scripts/common/`
`scripts/common/` holds code used by more than one data source (NOAA, gridMET,
NSRDB, Köppen). Scripts import from it, for example
`from common.counties import load_counties`; nothing in it is run directly.
Rules for `common/`:
- **Cross-source only.** A helper goes in only if metrics from more than one
data source use it. Code shared by scripts of a single data source stays with
that source, for example a future NSRDB module for the NSRDB prompt,
redaction, and error-log helpers.
- **One topic per module.** Each module is named for its topic and has a
docstring. No catch-all `utils.py`.
- **Keep it small.** Before adding a helper, ask why it does not belong to any
one data source.
Current modules:
| Module | Contents | Why it is in `common/` |
| --- | --- | --- |
| `county_zonal_stats.py` | Raster windows, the 180th-meridian split, area-weighted class shares | Used by any raster-based metric |
| `counties.py` | County polygon loading, FIPS normalization, the state FIPS table | County identity is shared by nearly every pipeline |
| `koppen_legend.py` | The Köppen code map and legend loader | Temporary: also used by `build_county_climate_data.py`; moves next to the Köppen code in Phase 3 |
**Rationale.** Helper functions make each step of a computation explicit,
avoid repeated code, and can be tested separately
([Brown CSCI 0111, "Helper Functions"](https://cs.brown.edu/courses/csci0111/fall2018/lectures/helper-functions.html)).
Shared helper folders, however, tend to lose cohesion and collect unrelated
code; the recommended alternative is to keep code with the part of the system
it belongs to, allowing a shared folder only if it stays small and documented
([Helpers and Utils Folders in Software Architecture](https://dev.to/knzt/helpers-and-utils-folders-in-software-architecture-3f8h)).
The rules above follow both: shared functions, organized by topic and limited
to code that crosses data sources.
Other shared code moves when its filter is reviewed, so each move is tested
alongside that filter. A 2026-09-13 survey found 19 functions with identical
copies in several scripts and 19 with copies that have drifted apart. Most
identical copies are NSRDB helpers, which belong in an NSRDB module rather
than `common/`; `read_csv_rows` (4 identical copies in `apply_*` scripts) is
cross-source. Drifted copies need a decision on which version is correct
before merging. Notable drifts: `summarize_county_gridmet_humidity.py` has its
own county loader and FIPS normalizer, and the state FIPS table is also copied
in `build_county_representative_points.py`,
`summarize_county_gridmet_humidity.py`, and
`request_nsrdb_county_polygon_archives.py`.
### Phase 3 — Assemble and restructure (after all 12 filters are clean)
- Assemble script that builds `climate-data.csv` from `data/metrics/`.
- Populate `metric_sources.json`; remove the per-row `source` column.
- Move audit columns out of the app CSV.
- Point the app's Sources panel at `metric_sources.json`.
- Retire or rewrite `build_county_climate_data.py` as per-metric builders.
- Move `common/koppen_legend.py` next to the Köppen code once nothing outside
Köppen imports it.
### Phase 4 — Runner and reproduction guide
- `scripts/pipeline.py` runner with `--only` and `--skip-download`.
- "Reproducing the data" guide organized by the tiers in section 3.
- Pinned requirements.
- End-to-end smoke test on a small synthetic county fixture.
- Organize scripts by data source (`noaa/`, `gridmet/`, `nsrdb/`, `koppen/`),
each holding its own helpers, with `common/` keeping only cross-source code.
Scripts in subfolders are run through the runner or as modules
(`python -m`), and the README and data-source commands are updated to match.
## 5. Per-filter review checklist
For each filter:
1. Verify the calculation against the source data and document findings.
2. Decide any rule or method changes with the project owner.
3. Record the adopted definition in `filter-calculations.md`.
4. Implement the calculation, writing `data/metrics/<metric>.csv`.
5. Add a single-column apply step for the current CSV.
6. Update the rules in `check_climate_data.py`.
7. Add or update unit tests.
8. Apply to the CSV, run `check_climate_data.py`, and compare changed counties
against expectations.
9. Update the app if the value set or display changes.
## 6. Filter tracker
Known issues come from `filter-calculations.md` ("Calculation review findings")
and this review; none beyond Köppen have been investigated yet.
| # | Filter | Status | Known issues to review |
| --- | --- | --- | --- |
| 1 | Köppen-Geiger class | Done (2026-09-14): rule applied, Mixed display built, documentation updated | See tasks below |
| 2 | Annual avg temperature | Not started | Months weighted equally (finding 1); touched-cell aggregation (finding 2) |
| 3 | Diurnal temperature range | Not started | Lexington, VA (51678) blank, while heat-index days use Rockbridge County as a proxy |
| 4 | Extreme temperature days | Not started | Partial years not normalized (finding 5); depends on retired percentile thresholds (finding 6); Lexington, VA blank |
| 5 | 90 °F+ heat-index days | Not started | Daily-extrema proxy (finding 7); permissive year completeness (finding 8) |
| 6 | Annual precipitation | Not started | Partial-year sums accepted (finding 3); touched-cell aggregation (finding 2) |
| 7 | Seasonality index | Not started | Touched-cell aggregation (finding 2) |
| 8 | Wettest month | Not started | Computed in both the base build and the precipitation-month script |
| 9 | Driest month | Not started | Same as wettest month |
| 10 | Summer specific humidity | Not started | Cell-center cos(latitude) aggregation differs from other metrics (finding 11) |
| 11 | Solar GHI | Not started | Finalized by the extreme-temperature apply script; hourly, 365-day assumption (finding 9) |
| 12 | Clear-sky GHI reduction | Not started | Mean of ratios rather than energy totals (finding 10) |
### Köppen-Geiger tasks
Adopted rule: a county is predominantly its top class if and only if that class
covers at least 50% of the county's land and leads the runner-up by at least
5 percentage points; otherwise it is Mixed climate. Expected result for the
50 states and DC: 3,010 predominant, 133 Mixed.
- [x] Investigate low-majority counties and adopt the rule.
- [x] Windowed raster reads and 180th-meridian split.
- [x] Area-weighted class shares (16 × 16 sub-cells per raster cell, scaled by
cos(latitude)) in `scripts/common/county_zonal_stats.py`.
- [x] Apply the 50% / 5-point rule (`scripts/build_county_koppen_metric.py`;
counties with no valid cells are left blank).
- [x] Run the builder to write `data/metrics/koppen.csv` and confirm the
expected 3,010 predominant / 133 Mixed (2026-09-13).
- [x] `koppenZone`-only apply step
(`scripts/apply_koppen_metric_to_climate_data.py`, with `--dry-run`).
- [x] Apply to `data/climate-data.csv` (2026-09-13; 142 counties changed to
Mixed, with `koppenPrimaryClass` and `koppenSecondaryClass` added).
- [x] Allow `Mixed` in `check_climate_data.py`.
- [x] Add a Mixed climate category to `app.js`, drawn as stripes of the
county's top two classes; see
[koppen-mixed-display-plan.md](koppen-mixed-display-plan.md).
- [x] Replace the plurality description in `filter-calculations.md` §1 and mark
review findings 2 and 4 resolved for Köppen (2026-09-14).
- [x] Update the Köppen descriptions and script lists in `README.md` and
`scripts/county_data_sources.md` (2026-09-14).
- [x] Tests for shares, the rule, boundary cases, and the apply step
(`tests/test_koppen_metric.py`).
## 7. Guardrails until Phase 3
- **Do not rerun `build_county_climate_data.py`** against
`data/climate-data.csv`. It would drop 9 live columns, restore
`extremeDays`, and overwrite the Mixed classification in `koppenZone`.
- Run `check_climate_data.py` after every apply step.
- Change one filter at a time, and compare its before and after values.
## 8. Open decisions
| Decision | Options | Needed by |
| --- | --- | --- |
| Committing large intermediates | Commit metric files only, or also source summaries | Phase 3 |
### Decided
- **Köppen audit columns (2026-09-12):** `koppen.csv` stores the top class and
share and the runner-up class and share alongside `koppenZone`.
- **Köppen no-data fallback (2026-09-12):** a county with no valid raster cells
is left blank, not assigned `Cfa`.
- **Applying Köppen to the app CSV (2026-09-12):** wait until the app supports
the Mixed class. Done 2026-09-13.
- **Metric files (2026-09-13):** one CSV per metric under `data/metrics/`,
starting with `koppen.csv`.
- **Mixed climate display (2026-09-13):** diagonal stripes of each Mixed
county's top two classes; see
[koppen-mixed-display-plan.md](koppen-mixed-display-plan.md).
- **Puerto Rico (2026-09-13):** off the map and out of every filter. The app
already drops state FIPS 72; the data files keep the rows.
- **Shared helpers (2026-09-13):** `scripts/common/` holds only code used by
more than one data source, one topic per module; code shared within one data
source stays with that source. Duplicates move during their own filter's
review.
+1 -1
View File
@@ -108,6 +108,6 @@
integrity="sha256-20nQCchB9co0qIjJZRGuk2/Z9VM+kNiyxNV1lvTlZBo=" integrity="sha256-20nQCchB9co0qIjJZRGuk2/Z9VM+kNiyxNV1lvTlZBo="
crossorigin="" crossorigin=""
></script> ></script>
<script src="app.js?v=precip-months-20260601"></script> <script src="app.js?v=koppen-mixed-cleanup-20260914"></script>
</body> </body>
</html> </html>
+15
View File
@@ -0,0 +1,15 @@
[tool.ruff]
line-length = 100
target-version = "py311"
[tool.ruff.lint]
select = ["B", "E", "F", "W", "I"]
ignore = ["E501"]
[tool.ruff.lint.per-file-ignores]
"tests/*.py" = ["E402"]
[tool.ruff.format]
quote-style = "double"
indent-style = "space"
line-ending = "auto"
@@ -14,7 +14,6 @@ from pathlib import Path
from build_county_locally_extreme_data import OLD_APP_FIPS_TO_CURRENT_FIPS from build_county_locally_extreme_data import OLD_APP_FIPS_TO_CURRENT_FIPS
REPO_ROOT = Path(__file__).resolve().parents[1] REPO_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv" DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv"
DEFAULT_DTR_CSV = REPO_ROOT / "data" / "noaa" / "county_diurnal_temperature_range.csv" DEFAULT_DTR_CSV = REPO_ROOT / "data" / "noaa" / "county_diurnal_temperature_range.csv"
@@ -15,7 +15,6 @@ import argparse
import csv import csv
from pathlib import Path from pathlib import Path
DEFAULT_CLIMATE_DATA = Path("data/climate-data.csv") DEFAULT_CLIMATE_DATA = Path("data/climate-data.csv")
DEFAULT_GRIDMET_HUMIDITY = Path("data/gridmet/county_gridmet_humidity.csv") DEFAULT_GRIDMET_HUMIDITY = Path("data/gridmet/county_gridmet_humidity.csv")
METRIC_FIELDS = [ METRIC_FIELDS = [
@@ -0,0 +1,138 @@
"""Apply the county Koppen-Geiger metric to the app CSV.
Replaces the koppenZone column of data/climate-data.csv with the values in
data/metrics/koppen.csv and writes koppenPrimaryClass and koppenSecondaryClass:
the two stripe classes the app draws for Mixed counties. Both are blank for
predominant counties. The two columns are added after koppenZone if missing.
Every other column, and the column order, is left unchanged. Use --dry-run to
report the changes without writing.
Run:
.venv\\Scripts\\python.exe scripts\\apply_koppen_metric_to_climate_data.py --dry-run
"""
from __future__ import annotations
import argparse
import csv
from collections import Counter
from pathlib import Path
from typing import Dict, List, Tuple
REPO_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv"
DEFAULT_KOPPEN_METRIC = REPO_ROOT / "data" / "metrics" / "koppen.csv"
ZONE_FIELD = "koppenZone"
PRIMARY_FIELD = "koppenPrimaryClass"
SECONDARY_FIELD = "koppenSecondaryClass"
STRIPE_FIELDS = [PRIMARY_FIELD, SECONDARY_FIELD]
MIXED_CLASS = "Mixed"
METRIC_FIELDS = ["countyFips", ZONE_FIELD, "koppenTopClass", "koppenSecondClass"]
# (countyFips, column, old value, new value)
Change = Tuple[str, str, str, str]
def read_csv_rows(path: Path) -> Tuple[List[str], List[dict]]:
"""Read a CSV while preserving the source field order."""
with path.open("r", encoding="utf-8-sig", newline="") as csv_file:
reader = csv.DictReader(csv_file)
if reader.fieldnames is None:
raise ValueError(f"{path} has no CSV header.")
return list(reader.fieldnames), list(reader)
def fieldnames_with_stripe_columns(fields: List[str]) -> List[str]:
"""Place the stripe-class columns directly after koppenZone."""
base = [field for field in fields if field not in STRIPE_FIELDS]
insert_at = base.index(ZONE_FIELD) + 1
return base[:insert_at] + STRIPE_FIELDS + base[insert_at:]
def load_koppen_values(koppen_metric: Path) -> Dict[str, Dict[str, str]]:
"""Return koppenZone and the two stripe classes for each county in the metric file."""
fields, rows = read_csv_rows(koppen_metric)
missing_fields = [field for field in METRIC_FIELDS if field not in fields]
if missing_fields:
raise ValueError(f"{koppen_metric} is missing columns {missing_fields}.")
values: Dict[str, Dict[str, str]] = {}
for row in rows:
is_mixed = row[ZONE_FIELD] == MIXED_CLASS
primary = row["koppenTopClass"] if is_mixed else ""
secondary = row["koppenSecondClass"] if is_mixed else ""
if is_mixed and not (primary and secondary):
raise ValueError(f"{row['countyFips']} is Mixed but has no top or second class in {koppen_metric}.")
values[row["countyFips"]] = {ZONE_FIELD: row[ZONE_FIELD], PRIMARY_FIELD: primary, SECONDARY_FIELD: secondary}
return values
def apply_koppen_metric(
climate_data: Path, koppen_metric: Path, out: Path, dry_run: bool = False
) -> Tuple[List[Change], List[str]]:
"""Update the Koppen columns; return the value changes and any columns that were added."""
fields, rows = read_csv_rows(climate_data)
if ZONE_FIELD not in fields:
raise ValueError(f"{climate_data} has no {ZONE_FIELD} column.")
lookup = load_koppen_values(koppen_metric)
missing = [row["countyFips"] for row in rows if row["countyFips"] not in lookup]
if missing:
raise ValueError(f"{len(missing)} counties are missing from {koppen_metric}, e.g. {missing[:5]}")
added_columns = [field for field in STRIPE_FIELDS if field not in fields]
changes: List[Change] = []
for row in rows:
new_values = lookup[row["countyFips"]]
for field in [ZONE_FIELD, *STRIPE_FIELDS]:
old_value = row.get(field) or ""
if old_value != new_values[field]:
changes.append((row["countyFips"], field, old_value, new_values[field]))
row[field] = new_values[field]
if not dry_run:
with out.open("w", encoding="utf-8", newline="") as csv_file:
writer = csv.DictWriter(csv_file, fieldnames=fieldnames_with_stripe_columns(fields))
writer.writeheader()
writer.writerows(rows)
return changes, added_columns
def parse_args() -> argparse.Namespace:
"""Define and parse command-line options for this apply step."""
parser = argparse.ArgumentParser(description="Apply the county Koppen metric to the app CSV.")
parser.add_argument("--climate-data", type=Path, default=DEFAULT_CLIMATE_DATA, help="App climate CSV to update.")
parser.add_argument("--koppen-metric", type=Path, default=DEFAULT_KOPPEN_METRIC, help="Koppen metric CSV.")
parser.add_argument("--out", type=Path, default=None, help="Output path; defaults to updating --climate-data in place.")
parser.add_argument("--dry-run", action="store_true", help="Report changes without writing.")
return parser.parse_args()
def main() -> None:
args = parse_args()
changes, added_columns = apply_koppen_metric(
args.climate_data, args.koppen_metric, args.out or args.climate_data, args.dry_run
)
verb = "would" if args.dry_run else "did"
if added_columns:
print(f"Columns added after {ZONE_FIELD} ({verb} write): {', '.join(added_columns)}")
zone_changes = [change for change in changes if change[1] == ZONE_FIELD]
kinds = Counter(
"to Mixed" if new == MIXED_CLASS else "to blank" if not new else "to another class"
for _, _, _, new in zone_changes
)
stripe_values = sum(1 for _, field, _, new in changes if field in STRIPE_FIELDS and new)
print(f"{len(zone_changes)} {ZONE_FIELD} values {'would change' if args.dry_run else 'changed'}: {dict(kinds)}")
print(f"{stripe_values} stripe-class values {'would be set' if args.dry_run else 'set'}.")
for fips, _, old, new in zone_changes[:10]:
print(f" {fips}: {old or '(blank)'} -> {new or '(blank)'}")
if len(zone_changes) > 10:
print(f" ... and {len(zone_changes) - 10} more")
if args.dry_run:
print("Dry run: nothing was written.")
if __name__ == "__main__":
main()
@@ -17,7 +17,6 @@ import argparse
import csv import csv
from pathlib import Path from pathlib import Path
REPO_ROOT = Path(__file__).resolve().parents[1] REPO_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv" DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv"
DEFAULT_COMPARISON = REPO_ROOT / "data" / "noaa" / "county_locally_extreme_days_comparison.csv" DEFAULT_COMPARISON = REPO_ROOT / "data" / "noaa" / "county_locally_extreme_days_comparison.csv"
@@ -84,13 +83,13 @@ def load_locally_extreme_lookup(comparison_csv: Path) -> dict[str, dict]:
def load_solar_ghi_lookup(solar_ghi_csv: Path) -> dict[str, str]: def load_solar_ghi_lookup(solar_ghi_csv: Path) -> dict[str, str]:
"""Map county FIPS codes to average daily GHI values from a county CSV.""" """Map county FIPS codes to mean daily GHI values from a county CSV."""
if not solar_ghi_csv.exists(): if not solar_ghi_csv.exists():
return {} return {}
fieldnames, rows = read_csv_rows(solar_ghi_csv) fieldnames, rows = read_csv_rows(solar_ghi_csv)
if "avgSolarGhiKwhM2Day" not in fieldnames: if "meanDailyGlobalHorizontalRadiationKwhM2Day" not in fieldnames:
raise ValueError(f"{solar_ghi_csv} must include avgSolarGhiKwhM2Day.") raise ValueError(f"{solar_ghi_csv} must include meanDailyGlobalHorizontalRadiationKwhM2Day.")
if "county_fips" in fieldnames: if "county_fips" in fieldnames:
fips_field = "county_fips" fips_field = "county_fips"
@@ -102,14 +101,14 @@ def load_solar_ghi_lookup(solar_ghi_csv: Path) -> dict[str, str]:
lookup: dict[str, str] = {} lookup: dict[str, str] = {}
for row in rows: for row in rows:
county_fips = normalize_fips(row.get(fips_field, "")) county_fips = normalize_fips(row.get(fips_field, ""))
raw_value = (row.get("avgSolarGhiKwhM2Day") or "").strip() raw_value = (row.get("meanDailyGlobalHorizontalRadiationKwhM2Day") or "").strip()
if not county_fips or not raw_value: if not county_fips or not raw_value:
continue continue
try: try:
float(raw_value) float(raw_value)
except ValueError as error: except ValueError as error:
raise ValueError( raise ValueError(
f"Invalid avgSolarGhiKwhM2Day value for county {county_fips}: {raw_value}" f"Invalid meanDailyGlobalHorizontalRadiationKwhM2Day value for county {county_fips}: {raw_value}"
) from error ) from error
lookup[county_fips] = raw_value lookup[county_fips] = raw_value
@@ -193,23 +192,23 @@ def apply_locally_extreme_metric(
polygon_solar_value = polygon_solar_lookup.get(county_fips) polygon_solar_value = polygon_solar_lookup.get(county_fips)
representative_solar_value = representative_solar_lookup.get(county_fips) representative_solar_value = representative_solar_lookup.get(county_fips)
current_solar_value = (row.get("avgSolarGhiKwhM2Day") or "").strip() current_solar_value = (row.get("meanDailyGlobalHorizontalRadiationKwhM2Day") or "").strip()
if polygon_solar_value: if polygon_solar_value:
row["avgSolarGhiKwhM2Day"] = polygon_solar_value row["meanDailyGlobalHorizontalRadiationKwhM2Day"] = polygon_solar_value
row["source"] = replace_solar_source_tag( row["source"] = replace_solar_source_tag(
row.get("source", ""), row.get("source", ""),
"solar-ghi-polygon-archive-area-weighted", "solar-ghi-polygon-archive-area-weighted",
) )
polygon_solar_count += 1 polygon_solar_count += 1
elif representative_solar_value or current_solar_value: elif representative_solar_value or current_solar_value:
row["avgSolarGhiKwhM2Day"] = representative_solar_value or current_solar_value row["meanDailyGlobalHorizontalRadiationKwhM2Day"] = representative_solar_value or current_solar_value
row["source"] = replace_solar_source_tag( row["source"] = replace_solar_source_tag(
row.get("source", ""), row.get("source", ""),
"solar-ghi-representative-point", "solar-ghi-representative-point",
) )
representative_solar_count += 1 representative_solar_count += 1
else: else:
row["avgSolarGhiKwhM2Day"] = "" row["meanDailyGlobalHorizontalRadiationKwhM2Day"] = ""
row["source"] = replace_solar_source_tag(row.get("source", ""), "missing-solar-ghi") row["source"] = replace_solar_source_tag(row.get("source", ""), "missing-solar-ghi")
missing_solar_count += 1 missing_solar_count += 1
@@ -221,9 +220,9 @@ def apply_locally_extreme_metric(
print(f"Wrote {len(climate_rows)} rows to {out}") print(f"Wrote {len(climate_rows)} rows to {out}")
print(f"Updated absoluteExtremeDays for {absolute_updated_count} rows.") print(f"Updated absoluteExtremeDays for {absolute_updated_count} rows.")
print(f"Left absoluteExtremeDays blank for {absolute_missing_count} rows without NOAA data.") print(f"Left absoluteExtremeDays blank for {absolute_missing_count} rows without NOAA data.")
print(f"Updated avgSolarGhiKwhM2Day from polygon archives for {polygon_solar_count} rows.") print(f"Updated meanDailyGlobalHorizontalRadiationKwhM2Day from polygon archives for {polygon_solar_count} rows.")
print(f"Kept representative-point solar fallback for {representative_solar_count} rows.") print(f"Kept representative-point solar fallback for {representative_solar_count} rows.")
print(f"Left avgSolarGhiKwhM2Day blank for {missing_solar_count} rows without solar data.") print(f"Left meanDailyGlobalHorizontalRadiationKwhM2Day blank for {missing_solar_count} rows without solar data.")
def parse_args() -> argparse.Namespace: def parse_args() -> argparse.Namespace:
@@ -1,9 +1,9 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
""" """
Merge NSRDB cloudiness metrics into climate-data.csv. Merge NSRDB clear-sky GHI reduction metrics into climate-data.csv.
Adds: Adds:
- cloudinessIndexPct - clearSkyGhiReductionIndex
Polygon area-weighted values are used first when available; representative-point Polygon area-weighted values are used first when available; representative-point
values remain the fallback. values remain the fallback.
@@ -15,33 +15,33 @@ import argparse
import csv import csv
from pathlib import Path from pathlib import Path
DEFAULT_CLIMATE_DATA = Path("data/climate-data.csv") DEFAULT_CLIMATE_DATA = Path("data/climate-data.csv")
DEFAULT_POLYGON_CLOUD_SUMMARY = Path("data/nrel/county_polygon_cloud_summary.csv") DEFAULT_POLYGON_CLOUD_SUMMARY = Path("data/nrel/county_polygon_cloud_summary.csv")
DEFAULT_REPRESENTATIVE_POINT_CLOUD_SUMMARY = Path("data/nrel/county_representative_point_cloud_summary.csv") DEFAULT_REPRESENTATIVE_POINT_CLOUD_SUMMARY = Path("data/nrel/county_representative_point_cloud_summary.csv")
METRIC_FIELD = "cloudinessIndexPct" METRIC_FIELD = "clearSkyGhiReductionIndex"
POLYGON_SOURCE_TAG = "nsrdb-polygon-area-weighted-cloudiness-tmy" POLYGON_METRIC_FIELD = "areaWeightedClearSkyGhiReductionIndex"
REPRESENTATIVE_POINT_SOURCE_TAG = "nsrdb-representative-point-cloudiness-tmy" POLYGON_SOURCE_TAG = "nsrdb-polygon-area-weighted-clear-sky-ghi-reduction-tmy"
REPRESENTATIVE_POINT_SOURCE_TAG = "nsrdb-representative-point-clear-sky-ghi-reduction-tmy"
CLOUD_SOURCE_TAGS = { CLOUD_SOURCE_TAGS = {
POLYGON_SOURCE_TAG, POLYGON_SOURCE_TAG,
REPRESENTATIVE_POINT_SOURCE_TAG, REPRESENTATIVE_POINT_SOURCE_TAG,
} }
def load_cloud_values(path: Path) -> dict[str, str]: def load_cloud_values(path: Path, metric_field: str = METRIC_FIELD) -> dict[str, str]:
"""Read cloudiness values keyed by county FIPS.""" """Read clear-sky GHI reduction values keyed by county FIPS."""
values: dict[str, str] = {} values: dict[str, str] = {}
with path.open("r", encoding="utf-8-sig", newline="") as handle: with path.open("r", encoding="utf-8-sig", newline="") as handle:
reader = csv.DictReader(handle) reader = csv.DictReader(handle)
fieldnames = reader.fieldnames or [] fieldnames = reader.fieldnames or []
required_fields = {"county_fips", METRIC_FIELD} required_fields = {"county_fips", metric_field}
missing_fields = sorted(required_fields - set(fieldnames)) missing_fields = sorted(required_fields - set(fieldnames))
if missing_fields: if missing_fields:
raise ValueError(f"{path} is missing fields: {', '.join(missing_fields)}") raise ValueError(f"{path} is missing fields: {', '.join(missing_fields)}")
for row in reader: for row in reader:
county_fips = (row.get("county_fips") or "").strip().zfill(5) county_fips = (row.get("county_fips") or "").strip().zfill(5)
value = (row.get(METRIC_FIELD) or "").strip() value = (row.get(metric_field) or "").strip()
if county_fips and value: if county_fips and value:
values[county_fips] = value values[county_fips] = value
return values return values
@@ -75,13 +75,13 @@ def merge_metric(
polygon_cloud_summary: Path, polygon_cloud_summary: Path,
representative_point_cloud_summary: Path, representative_point_cloud_summary: Path,
) -> tuple[int, int, int, int]: ) -> tuple[int, int, int, int]:
"""Merge cloudiness values into the app climate CSV.""" """Merge clear-sky GHI reduction values into the app climate CSV."""
polygon_cloudiness_by_fips = ( polygon_reduction_by_fips = (
load_cloud_values(polygon_cloud_summary) load_cloud_values(polygon_cloud_summary, POLYGON_METRIC_FIELD)
if polygon_cloud_summary.exists() if polygon_cloud_summary.exists()
else {} else {}
) )
representative_cloudiness_by_fips = load_cloud_values(representative_point_cloud_summary) representative_reduction_by_fips = load_cloud_values(representative_point_cloud_summary)
with climate_data.open("r", encoding="utf-8", newline="") as handle: with climate_data.open("r", encoding="utf-8", newline="") as handle:
reader = csv.DictReader(handle) reader = csv.DictReader(handle)
@@ -91,15 +91,15 @@ def merge_metric(
if "countyFips" not in fieldnames: if "countyFips" not in fieldnames:
raise ValueError(f"{climate_data} is missing countyFips") raise ValueError(f"{climate_data} is missing countyFips")
ensure_field_after(fieldnames, METRIC_FIELD, "avgSolarGhiKwhM2Day") ensure_field_after(fieldnames, METRIC_FIELD, "meanDailyGlobalHorizontalRadiationKwhM2Day")
polygon_count = 0 polygon_count = 0
representative_count = 0 representative_count = 0
missing_count = 0 missing_count = 0
for row in rows: for row in rows:
county_fips = (row.get("countyFips") or "").strip().zfill(5) county_fips = (row.get("countyFips") or "").strip().zfill(5)
polygon_value = polygon_cloudiness_by_fips.get(county_fips, "") polygon_value = polygon_reduction_by_fips.get(county_fips, "")
representative_value = representative_cloudiness_by_fips.get(county_fips, "") representative_value = representative_reduction_by_fips.get(county_fips, "")
value = polygon_value or representative_value value = polygon_value or representative_value
row[METRIC_FIELD] = value row[METRIC_FIELD] = value
if polygon_value: if polygon_value:
@@ -120,7 +120,7 @@ def merge_metric(
def main() -> int: def main() -> int:
parser = argparse.ArgumentParser(description="Merge NSRDB cloudiness metric into climate-data.csv.") parser = argparse.ArgumentParser(description="Merge NSRDB clear-sky GHI reduction metric into climate-data.csv.")
parser.add_argument("--climate-data", type=Path, default=DEFAULT_CLIMATE_DATA) parser.add_argument("--climate-data", type=Path, default=DEFAULT_CLIMATE_DATA)
parser.add_argument("--polygon-cloud-summary", type=Path, default=DEFAULT_POLYGON_CLOUD_SUMMARY) parser.add_argument("--polygon-cloud-summary", type=Path, default=DEFAULT_POLYGON_CLOUD_SUMMARY)
parser.add_argument( parser.add_argument(
@@ -14,16 +14,14 @@ from pathlib import Path
import numpy as np import numpy as np
import xarray as xr import xarray as xr
from build_county_climate_data import ( from build_county_climate_data import (
MONTH_NAMES, MONTH_NAMES,
_as_monthly_climatology, _as_monthly_climatology,
_extract_grid_2d, _extract_grid_2d,
_load_counties,
_select_data_var, _select_data_var,
_zonal_mean, _zonal_mean,
) )
from common.counties import load_counties
REPO_ROOT = Path(__file__).resolve().parents[1] REPO_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv" DEFAULT_CLIMATE_DATA = REPO_ROOT / "data" / "climate-data.csv"
@@ -63,7 +61,7 @@ def build_precip_month_lookup(
climatology_end_year: int, climatology_end_year: int,
) -> dict[str, tuple[str, str]]: ) -> dict[str, tuple[str, str]]:
"""Calculate wettest and driest precipitation month for each county.""" """Calculate wettest and driest precipitation month for each county."""
counties = _load_counties(counties_geojson) counties = load_counties(counties_geojson)
monthly_prcp = xr.open_dataset(monthly_prcp_nc, decode_times=True) monthly_prcp = xr.open_dataset(monthly_prcp_nc, decode_times=True)
try: try:
prcp_var = _select_data_var(monthly_prcp, "mlyprcp_norm") prcp_var = _select_data_var(monthly_prcp, "mlyprcp_norm")
+49 -220
View File
@@ -13,7 +13,7 @@ Metrics produced per county:
- wettestPrecipMonth: month with the highest 1991-2020 county mean precipitation - wettestPrecipMonth: month with the highest 1991-2020 county mean precipitation
- driestPrecipMonth: month with the lowest 1991-2020 county mean precipitation - driestPrecipMonth: month with the lowest 1991-2020 county mean precipitation
- extremeDays: count of normal-days with Tmax >= hot threshold or Tmin <= freeze threshold - extremeDays: count of normal-days with Tmax >= hot threshold or Tmin <= freeze threshold
- avgSolarGhiKwhM2Day: annual average daily global horizontal irradiance (GHI), when a solar raster or representative-point CSV is provided - meanDailyGlobalHorizontalRadiationKwhM2Day: mean daily global horizontal radiation (GHI), when a solar raster or representative-point CSV is provided
This script is intended for offline generation of complete county records. This script is intended for offline generation of complete county records.
""" """
@@ -22,111 +22,20 @@ from __future__ import annotations
import argparse import argparse
import csv import csv
import json
from pathlib import Path from pathlib import Path
from typing import Dict, List, Tuple from typing import Dict, List, Tuple
import geopandas as gpd import geopandas as gpd
import numpy as np import numpy as np
import rasterio import rasterio
from rasterio.features import geometry_mask
import xarray as xr import xarray as xr
from affine import Affine
from rasterio.features import geometry_mask
from shapely.geometry.base import BaseGeometry
DEFAULT_COUNTIES_GEOJSON_URL = "https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json" from common.counties import load_counties, normalize_fips
from common.county_zonal_stats import geometry_window, split_at_antimeridian
from common.koppen_legend import load_koppen_legend
STATE_FIPS_TO_ABBR = {
"01": "AL",
"02": "AK",
"04": "AZ",
"05": "AR",
"06": "CA",
"08": "CO",
"09": "CT",
"10": "DE",
"11": "DC",
"12": "FL",
"13": "GA",
"15": "HI",
"16": "ID",
"17": "IL",
"18": "IN",
"19": "IA",
"20": "KS",
"21": "KY",
"22": "LA",
"23": "ME",
"24": "MD",
"25": "MA",
"26": "MI",
"27": "MN",
"28": "MS",
"29": "MO",
"30": "MT",
"31": "NE",
"32": "NV",
"33": "NH",
"34": "NJ",
"35": "NM",
"36": "NY",
"37": "NC",
"38": "ND",
"39": "OH",
"40": "OK",
"41": "OR",
"42": "PA",
"44": "RI",
"45": "SC",
"46": "SD",
"47": "TN",
"48": "TX",
"49": "UT",
"50": "VT",
"51": "VA",
"53": "WA",
"54": "WV",
"55": "WI",
"56": "WY",
"60": "AS",
"66": "GU",
"69": "MP",
"72": "PR",
"78": "VI",
}
# Beck et al legend key is expected as text file, but this default handles common codes.
DEFAULT_KOPPEN_CODE_MAP = {
1: "Af",
2: "Am",
3: "Aw",
4: "BWh",
5: "BWk",
6: "BSh",
7: "BSk",
8: "Csa",
9: "Csb",
10: "Csc",
11: "Cwa",
12: "Cwb",
13: "Cwc",
14: "Cfa",
15: "Cfb",
16: "Cfc",
17: "Dsa",
18: "Dsb",
19: "Dsc",
20: "Dsd",
21: "Dwa",
22: "Dwb",
23: "Dwc",
24: "Dwd",
25: "Dfa",
26: "Dfb",
27: "Dfc",
28: "Dfd",
29: "ET",
30: "EF",
}
MONTH_NAMES = [ MONTH_NAMES = [
"January", "January",
@@ -144,94 +53,6 @@ MONTH_NAMES = [
] ]
def _normalize_fips(value: object, width: int) -> str:
"""Return a zero-padded FIPS code with the requested width."""
text = str(value).strip()
digits = "".join(ch for ch in text if ch.isdigit())
if not digits:
return ""
return digits.zfill(width)[-width:]
def _load_counties(counties_geojson: Path) -> gpd.GeoDataFrame:
"""Load county polygons and normalize fields used downstream."""
if not counties_geojson.exists():
try:
print(
f"County GeoJSON not found at {counties_geojson}. "
f"Attempting download from {DEFAULT_COUNTIES_GEOJSON_URL}..."
)
gdf = gpd.read_file(DEFAULT_COUNTIES_GEOJSON_URL)
counties_geojson.parent.mkdir(parents=True, exist_ok=True)
# Cache the downloaded file for subsequent runs.
gdf.to_file(counties_geojson, driver="GeoJSON")
print(f"Downloaded and cached county GeoJSON to {counties_geojson}")
except Exception as exc:
raise FileNotFoundError(
f"County GeoJSON not found at {counties_geojson}, and download from "
f"{DEFAULT_COUNTIES_GEOJSON_URL} failed. Download the file manually "
"and rerun with --counties-geojson pointing to it."
) from exc
gdf = gpd.read_file(counties_geojson)
if gdf.crs is None:
gdf = gdf.set_crs("EPSG:4326")
else:
gdf = gdf.to_crs("EPSG:4326")
feature_id = None
if "id" in gdf.columns:
feature_id = gdf["id"]
elif "GEOID" in gdf.columns:
feature_id = gdf["GEOID"]
elif "GEOID10" in gdf.columns:
feature_id = gdf["GEOID10"]
elif "fips" in gdf.columns:
feature_id = gdf["fips"]
else:
raise ValueError("Unable to locate county FIPS identifier column in county polygons.")
gdf["county_fips"] = feature_id.map(lambda value: _normalize_fips(value, 5))
gdf = gdf[gdf["county_fips"] != ""].copy()
if "NAME" in gdf.columns:
gdf["county_name"] = gdf["NAME"].fillna("").astype(str).str.strip()
elif "name" in gdf.columns:
gdf["county_name"] = gdf["name"].fillna("").astype(str).str.strip()
else:
gdf["county_name"] = gdf["county_fips"].map(lambda value: f"County {value}")
gdf["state_fips"] = gdf["county_fips"].str.slice(0, 2)
gdf["state"] = gdf["state_fips"].map(lambda code: STATE_FIPS_TO_ABBR.get(code, f"S{code}"))
gdf = gdf.sort_values("county_fips").reset_index(drop=True)
return gdf
def _load_koppen_legend(legend_path: Path | None) -> Dict[int, str]:
"""Load Koppen raster codes, using defaults when no legend exists."""
if legend_path is None:
return DEFAULT_KOPPEN_CODE_MAP
mapping: Dict[int, str] = {}
for line in legend_path.read_text(encoding="utf-8").splitlines():
text = line.strip()
if not text or text.startswith("#"):
continue
# Handles patterns like:
# "1: Af ..." or "1 = Af" or "1 Af"
import re
match = re.match(r"^(\d+)\s*[:=]?\s*([A-Za-z]{2,3})\b", text)
if not match:
continue
key = int(match.group(1))
value = match.group(2)
mapping[key] = value
return mapping if mapping else DEFAULT_KOPPEN_CODE_MAP
def _select_data_var(dataset: xr.Dataset, preferred: str) -> str: def _select_data_var(dataset: xr.Dataset, preferred: str) -> str:
"""Choose the best matching climate variable from a dataset.""" """Choose the best matching climate variable from a dataset."""
if preferred in dataset.data_vars: if preferred in dataset.data_vars:
@@ -261,7 +82,7 @@ def _select_data_var(dataset: xr.Dataset, preferred: str) -> str:
def _load_solar_ghi_csv(solar_ghi_csv: Path, counties: gpd.GeoDataFrame) -> List[float]: def _load_solar_ghi_csv(solar_ghi_csv: Path, counties: gpd.GeoDataFrame) -> List[float]:
"""Load county-keyed average daily GHI values from a representative-point or area-average CSV.""" """Load county-keyed mean daily GHI values from a representative-point or area-average CSV."""
with solar_ghi_csv.open(newline="", encoding="utf-8") as handle: with solar_ghi_csv.open(newline="", encoding="utf-8") as handle:
reader = csv.DictReader(handle) reader = csv.DictReader(handle)
if reader.fieldnames is None: if reader.fieldnames is None:
@@ -276,22 +97,22 @@ def _load_solar_ghi_csv(solar_ghi_csv: Path, counties: gpd.GeoDataFrame) -> List
f"Solar GHI CSV at {solar_ghi_csv} must include county_fips or countyFips." f"Solar GHI CSV at {solar_ghi_csv} must include county_fips or countyFips."
) )
if "avgSolarGhiKwhM2Day" not in reader.fieldnames: if "meanDailyGlobalHorizontalRadiationKwhM2Day" not in reader.fieldnames:
raise ValueError( raise ValueError(
f"Solar GHI CSV at {solar_ghi_csv} must include avgSolarGhiKwhM2Day." f"Solar GHI CSV at {solar_ghi_csv} must include meanDailyGlobalHorizontalRadiationKwhM2Day."
) )
solar_by_fips: Dict[str, float] = {} solar_by_fips: Dict[str, float] = {}
for row in reader: for row in reader:
county_fips = _normalize_fips(row.get(fips_field, ""), 5) county_fips = normalize_fips(row.get(fips_field, ""), 5)
raw_value = str(row.get("avgSolarGhiKwhM2Day", "")).strip() raw_value = str(row.get("meanDailyGlobalHorizontalRadiationKwhM2Day", "")).strip()
if not county_fips or not raw_value: if not county_fips or not raw_value:
continue continue
try: try:
solar_by_fips[county_fips] = float(raw_value) solar_by_fips[county_fips] = float(raw_value)
except ValueError as exc: except ValueError as exc:
raise ValueError( raise ValueError(
f"Invalid avgSolarGhiKwhM2Day value for county {county_fips}: {raw_value}" f"Invalid meanDailyGlobalHorizontalRadiationKwhM2Day value for county {county_fips}: {raw_value}"
) from exc ) from exc
return [ return [
@@ -357,7 +178,7 @@ def _infer_time_resolution_days(data_array: xr.DataArray) -> float:
return float(np.median(deltas)) return float(np.median(deltas))
def _extract_grid_2d(data_array: xr.DataArray) -> Tuple[np.ndarray, "Affine"]: def _extract_grid_2d(data_array: xr.DataArray) -> Tuple[np.ndarray, Affine]:
"""Convert a lat/lon slice to a raster array and transform.""" """Convert a lat/lon slice to a raster array and transform."""
# Expected shape for 2D arrays: lat, lon # Expected shape for 2D arrays: lat, lon
# Build affine from center coordinates. # Build affine from center coordinates.
@@ -376,8 +197,6 @@ def _extract_grid_2d(data_array: xr.DataArray) -> Tuple[np.ndarray, "Affine"]:
x_res = abs(lon[1] - lon[0]) x_res = abs(lon[1] - lon[0])
y_res = abs(lat[0] - lat[1]) y_res = abs(lat[0] - lat[1])
from affine import Affine
top_left_x = lon.min() - (x_res / 2.0) top_left_x = lon.min() - (x_res / 2.0)
top_left_y = lat.max() + (y_res / 2.0) top_left_y = lat.max() + (y_res / 2.0)
transform = Affine.translation(top_left_x, top_left_y) * Affine.scale(x_res, -y_res) transform = Affine.translation(top_left_x, top_left_y) * Affine.scale(x_res, -y_res)
@@ -416,6 +235,26 @@ def _zonal_mean_raster(raster_path: Path, counties: gpd.GeoDataFrame) -> List[fl
return _zonal_mean(values, source.transform, raster_counties) return _zonal_mean(values, source.transform, raster_counties)
def _touched_raster_values(source, geometry: BaseGeometry, split_antimeridian: bool) -> np.ndarray:
"""Read every raster cell a geometry touches, using a small window per piece."""
pieces = split_at_antimeridian(geometry) if split_antimeridian else [geometry]
selected: List[np.ndarray] = []
for piece in pieces:
window = geometry_window(source, piece)
values = source.read(1, window=window, masked=True).filled(0)
if source.nodata is not None:
values = np.where(values == source.nodata, 0, values)
mask = geometry_mask(
[piece.__geo_interface__],
out_shape=values.shape,
transform=source.window_transform(window),
invert=True,
all_touched=True,
)
selected.append(values[mask])
return np.concatenate(selected) if selected else np.array([], dtype=np.int64)
def _zonal_majority_class(koppen_raster: Path, counties: gpd.GeoDataFrame, code_map: Dict[int, str]) -> List[str]: def _zonal_majority_class(koppen_raster: Path, counties: gpd.GeoDataFrame, code_map: Dict[int, str]) -> List[str]:
"""Assign each county its most common Koppen-Geiger class.""" """Assign each county its most common Koppen-Geiger class."""
classes: List[str] = [] classes: List[str] = []
@@ -423,21 +262,11 @@ def _zonal_majority_class(koppen_raster: Path, counties: gpd.GeoDataFrame, code_
raster_counties = counties raster_counties = counties
if source.crs is not None and counties.crs is not None and counties.crs != source.crs: if source.crs is not None and counties.crs is not None and counties.crs != source.crs:
raster_counties = counties.to_crs(source.crs) raster_counties = counties.to_crs(source.crs)
# The 180-degree split only makes sense for longitude/latitude rasters.
data = source.read(1, masked=True) split_antimeridian = source.crs is None or source.crs.is_geographic
values = np.asarray(data.filled(0))
if source.nodata is not None:
values = np.where(values == source.nodata, 0, values)
for geometry in raster_counties.geometry: for geometry in raster_counties.geometry:
mask = geometry_mask( selected = _touched_raster_values(source, geometry, split_antimeridian)
[geometry.__geo_interface__],
out_shape=values.shape,
transform=source.transform,
invert=True,
all_touched=True,
)
selected = values[mask]
selected = selected[selected != 0] selected = selected[selected != 0]
if selected.size == 0: if selected.size == 0:
classes.append("Cfa") classes.append("Cfa")
@@ -477,7 +306,7 @@ def _compute_extreme_days(
day_tmax = _zonal_mean(tmax_arr, tmax_transform, counties) day_tmax = _zonal_mean(tmax_arr, tmax_transform, counties)
day_tmin = _zonal_mean(tmin_arr, tmin_transform, counties) day_tmin = _zonal_mean(tmin_arr, tmin_transform, counties)
for idx, (mx, mn) in enumerate(zip(day_tmax, day_tmin)): for idx, (mx, mn) in enumerate(zip(day_tmax, day_tmin, strict=True)):
if np.isnan(mx) or np.isnan(mn): if np.isnan(mx) or np.isnan(mn):
continue continue
if mx >= hot_threshold_c or mn <= freeze_threshold_c: if mx >= hot_threshold_c or mn <= freeze_threshold_c:
@@ -509,7 +338,7 @@ def _compute_extreme_days_monthly_proxy(
month_tmax = _zonal_mean(tmax_arr, tmax_transform, counties) month_tmax = _zonal_mean(tmax_arr, tmax_transform, counties)
month_tmin = _zonal_mean(tmin_arr, tmin_transform, counties) month_tmin = _zonal_mean(tmin_arr, tmin_transform, counties)
for idx, (mx, mn) in enumerate(zip(month_tmax, month_tmin)): for idx, (mx, mn) in enumerate(zip(month_tmax, month_tmin, strict=True)):
if np.isnan(mx) or np.isnan(mn): if np.isnan(mx) or np.isnan(mn):
continue continue
if mx >= hot_threshold_c or mn <= freeze_threshold_c: if mx >= hot_threshold_c or mn <= freeze_threshold_c:
@@ -547,9 +376,9 @@ def build_county_records(
solar_ghi_csv: Path | None, solar_ghi_csv: Path | None,
) -> Dict[str, dict]: ) -> Dict[str, dict]:
"""Build county climate records consumed by the web app.""" """Build county climate records consumed by the web app."""
counties = _load_counties(counties_geojson) counties = load_counties(counties_geojson)
koppen_classes = _zonal_majority_class(koppen_raster, counties, _load_koppen_legend(koppen_legend)) koppen_classes = _zonal_majority_class(koppen_raster, counties, load_koppen_legend(koppen_legend))
monthly_tavg = xr.open_dataset(monthly_tavg_nc, decode_times=True) monthly_tavg = xr.open_dataset(monthly_tavg_nc, decode_times=True)
monthly_prcp = xr.open_dataset(monthly_prcp_nc, decode_times=True) monthly_prcp = xr.open_dataset(monthly_prcp_nc, decode_times=True)
@@ -634,9 +463,9 @@ def build_county_records(
solar_source_tag = "solar-ghi-representative-point" solar_source_tag = "solar-ghi-representative-point"
else: else:
if solar_ghi_raster is not None: if solar_ghi_raster is not None:
print(f"Solar GHI raster not found at {solar_ghi_raster}; leaving avgSolarGhiKwhM2Day blank.") print(f"Solar GHI raster not found at {solar_ghi_raster}; leaving meanDailyGlobalHorizontalRadiationKwhM2Day blank.")
if solar_ghi_csv is not None: if solar_ghi_csv is not None:
print(f"Solar GHI CSV not found at {solar_ghi_csv}; leaving avgSolarGhiKwhM2Day blank.") print(f"Solar GHI CSV not found at {solar_ghi_csv}; leaving meanDailyGlobalHorizontalRadiationKwhM2Day blank.")
solar_ghi_kwh_m2_day = [float("nan")] * len(counties) solar_ghi_kwh_m2_day = [float("nan")] * len(counties)
solar_source_tag = "no-solar-ghi-source" solar_source_tag = "no-solar-ghi-source"
@@ -699,14 +528,14 @@ def build_county_records(
if not use_missing_extreme if not use_missing_extreme
else None else None
) )
avg_solar_ghi_value = ( mean_daily_global_horizontal_radiation_value = (
round(float(solar_ghi_kwh_m2_day[idx]), 2) round(float(solar_ghi_kwh_m2_day[idx]), 2)
if np.isfinite(solar_ghi_kwh_m2_day[idx]) if np.isfinite(solar_ghi_kwh_m2_day[idx])
else None else None
) )
source_suffix = " + missing-noaa-numeric" if used_any_missing_numeric else "" source_suffix = " + missing-noaa-numeric" if used_any_missing_numeric else ""
solar_source_suffix = "" if avg_solar_ghi_value is not None else " + missing-solar-ghi" solar_source_suffix = "" if mean_daily_global_horizontal_radiation_value is not None else " + missing-solar-ghi"
record = { record = {
"countyName": county_name, "countyName": county_name,
@@ -718,7 +547,7 @@ def build_county_records(
"wettestPrecipMonth": wettest_precip_month, "wettestPrecipMonth": wettest_precip_month,
"driestPrecipMonth": driest_precip_month, "driestPrecipMonth": driest_precip_month,
"extremeDays": extreme_days_value, "extremeDays": extreme_days_value,
"avgSolarGhiKwhM2Day": avg_solar_ghi_value, "meanDailyGlobalHorizontalRadiationKwhM2Day": mean_daily_global_horizontal_radiation_value,
"source": ( "source": (
"kg-beck2023 + noaa-nclimgrid-1991-2020 " "kg-beck2023 + noaa-nclimgrid-1991-2020 "
f"({extreme_days_source_tag}) + {solar_source_tag}{source_suffix}{solar_source_suffix}" f"({extreme_days_source_tag}) + {solar_source_tag}{source_suffix}{solar_source_suffix}"
@@ -751,7 +580,7 @@ def write_csv(records: Dict[str, dict], out_file: Path) -> None:
"wettestPrecipMonth", "wettestPrecipMonth",
"driestPrecipMonth", "driestPrecipMonth",
"extremeDays", "extremeDays",
"avgSolarGhiKwhM2Day", "meanDailyGlobalHorizontalRadiationKwhM2Day",
"source", "source",
] ]
@@ -818,14 +647,14 @@ def parse_args() -> argparse.Namespace:
"--solar-ghi-raster", "--solar-ghi-raster",
type=Path, type=Path,
default=None, default=None,
help="Optional raster of annual average daily GHI in kWh/m2/day for avgSolarGhiKwhM2Day.", help="Optional raster of mean daily GHI in kWh/m2/day for meanDailyGlobalHorizontalRadiationKwhM2Day.",
) )
parser.add_argument( parser.add_argument(
"--solar-ghi-csv", "--solar-ghi-csv",
type=Path, type=Path,
default=None, default=None,
help=( help=(
"Optional county CSV with avgSolarGhiKwhM2Day. Used as a representative-point fallback " "Optional county CSV with meanDailyGlobalHorizontalRadiationKwhM2Day. Used as a representative-point fallback "
"when --solar-ghi-raster is not supplied." "when --solar-ghi-raster is not supplied."
), ),
) )
@@ -27,7 +27,6 @@ from build_county_locally_extreme_data import (
rows_by_fips, rows_by_fips,
) )
REPO_ROOT = Path(__file__).resolve().parents[1] REPO_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_OUT = REPO_ROOT / "data" / "noaa" / "county_diurnal_temperature_range.csv" DEFAULT_OUT = REPO_ROOT / "data" / "noaa" / "county_diurnal_temperature_range.csv"
SOURCE_TAG = "noaa-nclimgrid-daily-county-area-averages-scaled" SOURCE_TAG = "noaa-nclimgrid-daily-county-area-averages-scaled"
@@ -90,7 +89,11 @@ def build_diurnal_temperature_range(
), ),
) )
for tmax_value, tmin_value in zip(tmax_record.values, tmin_record.values): for tmax_value, tmin_value in zip(
tmax_record.values,
tmin_record.values,
strict=True,
):
if tmax_value is None or tmin_value is None: if tmax_value is None or tmin_value is None:
continue continue
daily_range_c = tmax_value - tmin_value daily_range_c = tmax_value - tmin_value
+155
View File
@@ -0,0 +1,155 @@
"""Build the county Koppen-Geiger metric file from area-weighted class shares.
Writes data/metrics/koppen.csv. A county is predominantly its top class when
that class covers at least 50% of the county's land and leads the runner-up by
at least 5 percentage points; otherwise it is "Mixed". Counties with no valid
raster cells are left blank. See docs/filter-calculations.md, section 1.
Run:
.venv\\Scripts\\python.exe scripts\\build_county_koppen_metric.py
"""
from __future__ import annotations
import argparse
import csv
from collections import Counter
from pathlib import Path
from typing import Dict, List, Tuple
import geopandas as gpd
import rasterio
from common.counties import load_counties
from common.county_zonal_stats import DEFAULT_SUBCELLS, area_weighted_class_weights
from common.koppen_legend import load_koppen_legend
PROJECT_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_COUNTIES_GEOJSON = PROJECT_ROOT / "data" / "geojson-counties-fips.json"
DEFAULT_KOPPEN_RASTER = PROJECT_ROOT / "data" / "koppen_geiger_tif" / "1991_2020" / "koppen_geiger_0p00833333.tif"
DEFAULT_KOPPEN_LEGEND = PROJECT_ROOT / "data" / "koppen_geiger_tif" / "legend.txt"
DEFAULT_OUT = PROJECT_ROOT / "data" / "metrics" / "koppen.csv"
PREDOMINANT_MIN_SHARE = 0.50
PREDOMINANT_MIN_GAP = 0.05
MIXED_CLASS = "Mixed"
# Keeps shares that sit exactly on a cutoff from failing on floating-point error.
RULE_TOLERANCE = 1e-9
FIELDS = [
"countyFips",
"countyName",
"state",
"koppenZone",
"koppenTopClass",
"koppenTopShare",
"koppenSecondClass",
"koppenSecondShare",
]
def rank_class_shares(weights: Dict[int, float], code_map: Dict[int, str]) -> List[Tuple[str, float]]:
"""Convert per-code area weights to class shares, largest first.
Exact ties go to the smaller raster code, matching the previous build.
"""
total = sum(weights.values())
if total <= 0:
return []
unknown = sorted(code for code in weights if code not in code_map)
if unknown:
raise ValueError(f"Raster codes {unknown} are not in the Koppen legend.")
ranked = sorted(weights.items(), key=lambda item: (-item[1], item[0]))
return [(code_map[code], weight / total) for code, weight in ranked]
def classify(ranked: List[Tuple[str, float]]) -> str:
"""Return the predominant class, "Mixed", or blank when there is no data."""
if not ranked:
return ""
top_share = ranked[0][1]
second_share = ranked[1][1] if len(ranked) > 1 else 0.0
has_majority = top_share >= PREDOMINANT_MIN_SHARE - RULE_TOLERANCE
has_clear_lead = top_share - second_share >= PREDOMINANT_MIN_GAP - RULE_TOLERANCE
return ranked[0][0] if has_majority and has_clear_lead else MIXED_CLASS
def _format_share(ranked: List[Tuple[str, float]], index: int) -> Tuple[str, str]:
"""Return the class and 4-decimal share at a rank, or blanks."""
if index >= len(ranked):
return "", ""
code, share = ranked[index]
return code, f"{share:.4f}"
def build_koppen_records(
counties: gpd.GeoDataFrame,
koppen_raster: Path,
code_map: Dict[int, str],
subcells: int = DEFAULT_SUBCELLS,
) -> List[dict]:
"""Classify every county and return its metric row."""
records: List[dict] = []
with rasterio.open(koppen_raster) as source:
raster_counties = counties
if source.crs is not None and counties.crs is not None and counties.crs != source.crs:
raster_counties = counties.to_crs(source.crs)
geographic = source.crs is None or source.crs.is_geographic
for county, geometry in zip(counties.itertuples(), raster_counties.geometry):
weights = area_weighted_class_weights(source, geometry, geographic=geographic, subcells=subcells)
ranked = rank_class_shares(weights, code_map)
top_class, top_share = _format_share(ranked, 0)
second_class, second_share = _format_share(ranked, 1)
records.append(
{
"countyFips": county.county_fips,
"countyName": county.county_name,
"state": county.state,
"koppenZone": classify(ranked),
"koppenTopClass": top_class,
"koppenTopShare": top_share,
"koppenSecondClass": second_class,
"koppenSecondShare": second_share,
}
)
return records
def write_records(records: List[dict], out_file: Path) -> None:
"""Write the Koppen metric rows."""
out_file.parent.mkdir(parents=True, exist_ok=True)
with out_file.open("w", encoding="utf-8", newline="") as csv_file:
writer = csv.DictWriter(csv_file, fieldnames=FIELDS)
writer.writeheader()
writer.writerows(records)
def parse_args() -> argparse.Namespace:
"""Define and parse command-line options for this builder."""
parser = argparse.ArgumentParser(description="Build the county Koppen-Geiger metric file.")
parser.add_argument("--counties-geojson", type=Path, default=DEFAULT_COUNTIES_GEOJSON, help="County polygon GeoJSON path.")
parser.add_argument("--koppen-raster", type=Path, default=DEFAULT_KOPPEN_RASTER, help="Koppen-Geiger raster TIFF path.")
parser.add_argument("--koppen-legend", type=Path, default=DEFAULT_KOPPEN_LEGEND, help="legend.txt mapping raster codes.")
parser.add_argument("--subcells", type=int, default=DEFAULT_SUBCELLS, help="Sub-cells per raster cell edge.")
parser.add_argument("--out", type=Path, default=DEFAULT_OUT, help="Output metric CSV path.")
return parser.parse_args()
def main() -> None:
args = parse_args()
counties = load_counties(args.counties_geojson)
records = build_koppen_records(counties, args.koppen_raster, load_koppen_legend(args.koppen_legend), args.subcells)
write_records(records, args.out)
outcomes = Counter(
"blank" if not r["koppenZone"] else "Mixed" if r["koppenZone"] == MIXED_CLASS else "predominant" for r in records
)
print(
f"Wrote {len(records)} counties to {args.out}: {outcomes['predominant']} predominant, "
f"{outcomes['Mixed']} Mixed, {outcomes['blank']} blank."
)
if __name__ == "__main__":
main()
+5 -2
View File
@@ -36,7 +36,6 @@ from pathlib import Path
from tempfile import NamedTemporaryFile from tempfile import NamedTemporaryFile
from typing import Dict, Iterable, Iterator, List, Tuple from typing import Dict, Iterable, Iterator, List, Tuple
NOAA_AVERAGES_BASE_URL = "https://www.ncei.noaa.gov/data/nclimgrid-daily/access/averages" NOAA_AVERAGES_BASE_URL = "https://www.ncei.noaa.gov/data/nclimgrid-daily/access/averages"
NOAA_STATE_CROSSWALK_URL = ( NOAA_STATE_CROSSWALK_URL = (
"https://www.ncei.noaa.gov/data/nclimgrid-daily/doc/us-state-codes_ncei-to-fips.csv" "https://www.ncei.noaa.gov/data/nclimgrid-daily/doc/us-state-codes_ncei-to-fips.csv"
@@ -518,7 +517,11 @@ def build_annual_counts(
) )
counts = annual_counts.setdefault((county_fips, year), AnnualCounts()) counts = annual_counts.setdefault((county_fips, year), AnnualCounts())
for tmax_value, tmin_value in zip(tmax_record.values, tmin_record.values): for tmax_value, tmin_value in zip(
tmax_record.values,
tmin_record.values,
strict=True,
):
if tmax_value is None or tmin_value is None: if tmax_value is None or tmin_value is None:
continue continue
@@ -13,7 +13,6 @@ from pathlib import Path
import geopandas as gpd import geopandas as gpd
DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json") DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json")
DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_points.csv") DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_points.csv")
+427
View File
@@ -0,0 +1,427 @@
"""Validate data/climate-data.csv before the browser app loads it.
Checks file format, columns, county keys, the county geometry join, per-metric
types and ranges, allowed blank values, and cross-field consistency. Exits with
status 1 when any check fails.
Run:
.venv\\Scripts\\python.exe scripts\\check_climate_data.py
"""
from __future__ import annotations
import argparse
import csv
import json
import math
import re
import sys
from dataclasses import dataclass, field
from pathlib import Path
from typing import Dict, FrozenSet, List, Optional, Sequence, Tuple
PROJECT_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_CLIMATE_CSV = PROJECT_ROOT / "data" / "climate-data.csv"
DEFAULT_COUNTIES_GEOJSON = PROJECT_ROOT / "data" / "geojson-counties-fips.json"
DEFAULT_METRIC_SOURCES = PROJECT_ROOT / "data" / "metric_sources.json"
FIPS_PATTERN = re.compile(r"\d{5}")
STATE_PATTERN = re.compile(r"[A-Z]{2}")
# The distinct-value check only runs on full-size files, not small fixtures.
COLLAPSE_CHECK_MIN_ROWS = 100
KOPPEN_CODES = (
"Af", "Am", "Aw", "BWh", "BWk", "BSh", "BSk",
"Csa", "Csb", "Csc", "Cwa", "Cwb", "Cwc", "Cfa", "Cfb", "Cfc",
"Dsa", "Dsb", "Dsc", "Dsd", "Dwa", "Dwb", "Dwc", "Dwd",
"Dfa", "Dfb", "Dfc", "Dfd", "ET", "EF",
)
# Counties with no predominant class; see docs/filter-calculations.md, section 1.
MIXED_KOPPEN_CLASS = "Mixed"
MONTH_NAMES = (
"January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December",
)
# NOAA nClimGrid and gridMET cover the contiguous U.S. only.
OUTSIDE_CONUS = frozenset({"AK", "HI", "PR"})
# NSRDB summaries were not requested for Puerto Rico.
NSRDB_NOT_REQUESTED = frozenset({"PR"})
# NOAA county daily files do not include Lexington city, VA.
NOAA_DAILY_MISSING_FIPS = frozenset({"51678"})
@dataclass(frozen=True)
class MetricRule:
"""Validation rule for one metric column."""
kind: str
minimum: Optional[float] = None
maximum: Optional[float] = None
integer: bool = False
categories: Tuple[str, ...] = ()
blank_states: FrozenSet[str] = frozenset()
blank_fips: FrozenSet[str] = frozenset()
min_distinct: int = 2
def _numeric(
minimum: float,
maximum: float,
*,
integer: bool = False,
blank_states: FrozenSet[str] = frozenset(),
blank_fips: FrozenSet[str] = frozenset(),
) -> MetricRule:
"""Build a numeric rule with a plausible physical range."""
return MetricRule(
"numeric",
minimum=minimum,
maximum=maximum,
integer=integer,
blank_states=blank_states,
blank_fips=blank_fips,
min_distinct=10,
)
def _categorical(categories: Sequence[str], *, blank_states: FrozenSet[str] = frozenset()) -> MetricRule:
"""Build a categorical rule limited to an allowed value list."""
return MetricRule("categorical", categories=tuple(categories), blank_states=blank_states)
# Ranges are physical plausibility limits, deliberately wider than the current
# data. They are not the app's color-scale bounds.
METRIC_RULES: Dict[str, MetricRule] = {
"koppenZone": _categorical(KOPPEN_CODES + (MIXED_KOPPEN_CLASS,)),
"avgTempF": _numeric(20, 85, blank_states=OUTSIDE_CONUS),
"avgDiurnalTempRangeF": _numeric(5, 40, blank_states=OUTSIDE_CONUS, blank_fips=NOAA_DAILY_MISSING_FIPS),
"annualPrecipIn": _numeric(1, 150, blank_states=OUTSIDE_CONUS),
"seasonalityIndex": _numeric(0, 100, integer=True, blank_states=OUTSIDE_CONUS),
"wettestPrecipMonth": _categorical(MONTH_NAMES, blank_states=OUTSIDE_CONUS),
"driestPrecipMonth": _categorical(MONTH_NAMES, blank_states=OUTSIDE_CONUS),
"absoluteExtremeDays": _numeric(0, 366, blank_states=OUTSIDE_CONUS, blank_fips=NOAA_DAILY_MISSING_FIPS),
"meanDailyGlobalHorizontalRadiationKwhM2Day": _numeric(1.5, 7.5, blank_states=NSRDB_NOT_REQUESTED),
"clearSkyGhiReductionIndex": _numeric(0, 1, blank_states=NSRDB_NOT_REQUESTED),
"avgSummerSpecificHumidityGKg": _numeric(1, 25, blank_states=OUTSIDE_CONUS),
"humidHeatDays": _numeric(0, 366, blank_states=OUTSIDE_CONUS),
}
IDENTITY_COLUMNS = ("countyFips", "countyName", "state")
# Stripe classes the app draws for Mixed Koppen counties; blank otherwise.
KOPPEN_STRIPE_COLUMNS = ("koppenPrimaryClass", "koppenSecondaryClass")
AUDIT_COLUMNS = ("humidHeatSourceFips", "humidHeatFipsAdjustment", "source")
EXPECTED_COLUMNS = IDENTITY_COLUMNS + tuple(METRIC_RULES) + KOPPEN_STRIPE_COLUMNS + AUDIT_COLUMNS
CHECK_FORMAT = "file format"
CHECK_COLUMNS = "columns"
CHECK_COUNTIES = "county keys"
CHECK_GEOMETRY = "map geometry join"
CHECK_BLANKS = "blank values"
CHECK_VALUES = "value types and ranges"
CHECK_CROSS = "cross-field consistency"
CHECK_SOURCES = "metric sources file"
@dataclass
class Report:
"""Problems found per check, in the order checks ran."""
checks: Dict[str, List[str]] = field(default_factory=dict)
notes: List[str] = field(default_factory=list)
row_count: int = 0
column_count: int = 0
def start(self, check: str) -> None:
self.checks.setdefault(check, [])
def add(self, check: str, message: str) -> None:
self.checks.setdefault(check, []).append(message)
@property
def problem_count(self) -> int:
return sum(len(messages) for messages in self.checks.values())
@property
def ok(self) -> bool:
return self.problem_count == 0
def read_climate_csv(path: Path) -> Tuple[List[str], List[List[str]]]:
"""Read the raw header and rows without DictReader's silent padding."""
with path.open("r", encoding="utf-8", newline="") as handle:
reader = csv.reader(handle)
headers = next(reader, [])
rows = [row for row in reader if row]
return headers, rows
def check_structure(
report: Report, headers: List[str], raw_rows: List[List[str]]
) -> Tuple[List[str], List[Dict[str, str]]]:
"""Check encoding, header names, and row widths; return rows as dicts."""
report.start(CHECK_FORMAT)
report.start(CHECK_COLUMNS)
if headers and headers[0].startswith("\ufeff"):
report.add(CHECK_FORMAT, "file starts with a UTF-8 BOM; the app would not find the first column")
headers = [headers[0].lstrip("\ufeff")] + headers[1:]
seen = set()
for name in headers:
if name in seen:
report.add(CHECK_COLUMNS, f"duplicate column {name!r}")
seen.add(name)
for name in EXPECTED_COLUMNS:
if name not in seen:
report.add(CHECK_COLUMNS, f"missing column {name!r}")
for name in headers:
if name not in EXPECTED_COLUMNS:
report.add(CHECK_COLUMNS, f"unexpected column {name!r} (add it to check_climate_data.py if intended)")
rows: List[Dict[str, str]] = []
for index, raw in enumerate(raw_rows):
if len(raw) != len(headers):
report.add(CHECK_FORMAT, f"line {index + 2}: {len(raw)} fields, expected {len(headers)}")
continue
rows.append(dict(zip(headers, raw)))
report.row_count = len(raw_rows)
report.column_count = len(headers)
return headers, rows
def check_counties(report: Report, rows: List[Dict[str, str]]) -> None:
"""Check FIPS format and uniqueness, names, and state/FIPS agreement."""
report.start(CHECK_COUNTIES)
seen = set()
prefix_states: Dict[str, set] = {}
state_prefixes: Dict[str, set] = {}
for row in rows:
fips = row.get("countyFips") or ""
if not FIPS_PATTERN.fullmatch(fips):
report.add(CHECK_COUNTIES, f"invalid countyFips {fips!r}")
continue
if fips in seen:
report.add(CHECK_COUNTIES, f"duplicate countyFips {fips}")
seen.add(fips)
if not (row.get("countyName") or "").strip():
report.add(CHECK_COUNTIES, f"{fips}: countyName is blank")
state = row.get("state") or ""
if not STATE_PATTERN.fullmatch(state):
report.add(CHECK_COUNTIES, f"{fips}: invalid state {state!r}")
continue
prefix_states.setdefault(fips[:2], set()).add(state)
state_prefixes.setdefault(state, set()).add(fips[:2])
for prefix, states in sorted(prefix_states.items()):
if len(states) > 1:
report.add(CHECK_COUNTIES, f"state FIPS {prefix} is labeled as {sorted(states)}")
for state, prefixes in sorted(state_prefixes.items()):
if len(prefixes) > 1:
report.add(CHECK_COUNTIES, f"state {state} spans FIPS prefixes {sorted(prefixes)}")
def _feature_fips(feature: dict) -> str:
"""Return the 5-digit county FIPS stored on a GeoJSON feature."""
props = feature.get("properties") or {}
raw = feature.get("id") or props.get("id") or str(props.get("GEO_ID") or "")[-5:]
text = str(raw or "").strip()
return text.zfill(5) if text.isdigit() else text
def check_geometry(report: Report, rows: List[Dict[str, str]], geojson_path: Path) -> None:
"""Check that CSV counties and map polygons match one-to-one."""
report.start(CHECK_GEOMETRY)
if not geojson_path.exists():
report.add(CHECK_GEOMETRY, f"county geometry file not found: {geojson_path}")
return
features = json.loads(geojson_path.read_text(encoding="utf-8")).get("features", [])
geo_fips = set()
for feature in features:
fips = _feature_fips(feature)
if not FIPS_PATTERN.fullmatch(fips):
report.add(CHECK_GEOMETRY, f"map feature without a valid county FIPS: {fips!r}")
continue
if fips in geo_fips:
report.add(CHECK_GEOMETRY, f"map has duplicate polygons for {fips}")
geo_fips.add(fips)
state_fips = str((feature.get("properties") or {}).get("STATE") or "").strip()
if state_fips and state_fips.zfill(2) != fips[:2]:
report.add(CHECK_GEOMETRY, f"map feature {fips} has STATE {state_fips}")
csv_fips = {row.get("countyFips") or "" for row in rows}
for fips in sorted(csv_fips - geo_fips):
report.add(CHECK_GEOMETRY, f"{fips} is in the CSV but has no map polygon")
for fips in sorted(geo_fips - csv_fips):
report.add(CHECK_GEOMETRY, f"{fips} has a map polygon but no CSV row")
def check_metrics(report: Report, headers: List[str], rows: List[Dict[str, str]]) -> None:
"""Check each metric's blanks, types, ranges, and categories."""
report.start(CHECK_BLANKS)
report.start(CHECK_VALUES)
for key, rule in METRIC_RULES.items():
if key not in headers:
continue
present: List[object] = []
for row in rows:
fips = row.get("countyFips") or "?"
state = row.get("state") or ""
value = row.get(key) or ""
if value == "":
if state not in rule.blank_states and fips not in rule.blank_fips:
report.add(CHECK_BLANKS, f"{fips} ({state}): {key} is blank")
continue
if value != value.strip():
report.add(CHECK_VALUES, f"{fips}: {key}={value!r} has surrounding whitespace")
continue
if rule.kind == "categorical":
if value in rule.categories:
present.append(value)
else:
report.add(CHECK_VALUES, f"{fips}: {key}={value!r} is not an allowed category")
continue
try:
number = float(value)
except ValueError:
report.add(CHECK_VALUES, f"{fips}: {key}={value!r} is not a number")
continue
if not math.isfinite(number):
report.add(CHECK_VALUES, f"{fips}: {key}={value!r} is not finite")
continue
if rule.integer and not number.is_integer():
report.add(CHECK_VALUES, f"{fips}: {key}={value} should be a whole number")
if not rule.minimum <= number <= rule.maximum:
report.add(CHECK_VALUES, f"{fips}: {key}={value} outside {rule.minimum:g}..{rule.maximum:g}")
present.append(number)
if len(rows) >= COLLAPSE_CHECK_MIN_ROWS and len(set(present)) < rule.min_distinct:
report.add(
CHECK_VALUES,
f"{key} has only {len(set(present))} distinct values; the column may have been overwritten",
)
def check_cross_fields(report: Report, rows: List[Dict[str, str]]) -> None:
"""Check relationships between columns within each row."""
report.start(CHECK_CROSS)
for row in rows:
fips = row.get("countyFips") or "?"
wet = row.get("wettestPrecipMonth") or ""
dry = row.get("driestPrecipMonth") or ""
if bool(wet) != bool(dry):
report.add(CHECK_CROSS, f"{fips}: only one of wettestPrecipMonth/driestPrecipMonth is set")
elif wet and wet == dry:
report.add(CHECK_CROSS, f"{fips}: wettest and driest month are both {wet}")
heat = row.get("humidHeatDays") or ""
heat_source = row.get("humidHeatSourceFips") or ""
if bool(heat) != bool(heat_source):
report.add(CHECK_CROSS, f"{fips}: humidHeatDays and humidHeatSourceFips must both be set or both blank")
if heat_source and not FIPS_PATTERN.fullmatch(heat_source):
report.add(CHECK_CROSS, f"{fips}: invalid humidHeatSourceFips {heat_source!r}")
elif heat_source and heat_source != fips and not (row.get("humidHeatFipsAdjustment") or "").strip():
report.add(CHECK_CROSS, f"{fips}: uses proxy county {heat_source} without a humidHeatFipsAdjustment note")
zone = row.get("koppenZone") or ""
primary = row.get("koppenPrimaryClass") or ""
secondary = row.get("koppenSecondaryClass") or ""
if zone == MIXED_KOPPEN_CLASS:
if not (primary and secondary):
report.add(CHECK_CROSS, f"{fips}: Mixed Koppen county needs koppenPrimaryClass and koppenSecondaryClass")
elif primary not in KOPPEN_CODES or secondary not in KOPPEN_CODES:
report.add(CHECK_CROSS, f"{fips}: invalid Koppen stripe classes {primary!r}/{secondary!r}")
elif primary == secondary:
report.add(CHECK_CROSS, f"{fips}: Koppen stripe classes are both {primary}")
elif primary or secondary:
report.add(CHECK_CROSS, f"{fips}: Koppen stripe classes are set but koppenZone is {zone or 'blank'}, not Mixed")
if not (row.get("source") or "").strip():
report.add(CHECK_CROSS, f"{fips}: source is blank")
def check_metric_sources(report: Report, path: Path) -> None:
"""Check that the metric metadata file parses and names real metrics."""
report.start(CHECK_SOURCES)
if not path.exists():
report.add(CHECK_SOURCES, f"{path.name} not found")
return
try:
data = json.loads(path.read_text(encoding="utf-8"))
except json.JSONDecodeError as exc:
report.add(CHECK_SOURCES, f"{path.name} is not valid JSON: {exc}")
return
if not isinstance(data, dict) or not isinstance(data.get("metrics"), dict):
report.add(CHECK_SOURCES, f'{path.name} must be a JSON object with a "metrics" object')
return
documented = set(data["metrics"])
for key in sorted(documented - set(METRIC_RULES)):
report.add(CHECK_SOURCES, f"{path.name} describes unknown metric {key!r}")
covered = len(documented & set(METRIC_RULES))
report.notes.append(f"{path.name} documents {covered} of {len(METRIC_RULES)} metrics.")
def run_checks(csv_path: Path, geojson_path: Path, metric_sources_path: Path) -> Report:
"""Run every check and return the combined report."""
report = Report()
headers, raw_rows = read_climate_csv(csv_path)
headers, rows = check_structure(report, headers, raw_rows)
check_counties(report, rows)
check_geometry(report, rows, geojson_path)
check_metrics(report, headers, rows)
check_cross_fields(report, rows)
check_metric_sources(report, metric_sources_path)
return report
def print_report(report: Report, csv_path: Path, max_examples: int) -> None:
"""Print a pass/fail line per check with example problems."""
print(f"Checked {csv_path}: {report.row_count} rows, {report.column_count} columns")
for check, messages in report.checks.items():
status = "PASS" if not messages else f"FAIL ({len(messages)})"
print(f" {status:<11}{check}")
for message in messages[:max_examples]:
print(f" - {message}")
if len(messages) > max_examples:
print(f" ... and {len(messages) - max_examples} more")
for note in report.notes:
print(f" note: {note}")
print("Result: PASS" if report.ok else f"Result: FAIL ({report.problem_count} problems)")
def parse_args() -> argparse.Namespace:
"""Define and parse command-line options for this checker."""
parser = argparse.ArgumentParser(description="Validate the browser app's county climate CSV.")
parser.add_argument("--csv", type=Path, default=DEFAULT_CLIMATE_CSV, help="Climate data CSV path.")
parser.add_argument("--geojson", type=Path, default=DEFAULT_COUNTIES_GEOJSON, help="County GeoJSON path.")
parser.add_argument("--metric-sources", type=Path, default=DEFAULT_METRIC_SOURCES, help="Metric metadata JSON path.")
parser.add_argument("--max-examples", type=int, default=10, help="Problems to print per failing check.")
return parser.parse_args()
def main() -> int:
args = parse_args()
if not args.csv.exists():
print(f"Climate data CSV not found: {args.csv}", file=sys.stderr)
return 2
report = run_checks(args.csv, args.geojson, args.metric_sources)
print_report(report, args.csv, args.max_examples)
return 0 if report.ok else 1
if __name__ == "__main__":
sys.exit(main())
+5
View File
@@ -0,0 +1,5 @@
"""Shared helpers imported by the county data pipeline scripts.
Modules here are not run directly. Scripts in the parent folder import them,
for example ``from common.counties import load_counties``.
"""
+131
View File
@@ -0,0 +1,131 @@
"""Load county polygons and normalize county identifiers."""
from __future__ import annotations
from pathlib import Path
import geopandas as gpd
DEFAULT_COUNTIES_GEOJSON_URL = "https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json"
STATE_FIPS_TO_ABBR = {
"01": "AL",
"02": "AK",
"04": "AZ",
"05": "AR",
"06": "CA",
"08": "CO",
"09": "CT",
"10": "DE",
"11": "DC",
"12": "FL",
"13": "GA",
"15": "HI",
"16": "ID",
"17": "IL",
"18": "IN",
"19": "IA",
"20": "KS",
"21": "KY",
"22": "LA",
"23": "ME",
"24": "MD",
"25": "MA",
"26": "MI",
"27": "MN",
"28": "MS",
"29": "MO",
"30": "MT",
"31": "NE",
"32": "NV",
"33": "NH",
"34": "NJ",
"35": "NM",
"36": "NY",
"37": "NC",
"38": "ND",
"39": "OH",
"40": "OK",
"41": "OR",
"42": "PA",
"44": "RI",
"45": "SC",
"46": "SD",
"47": "TN",
"48": "TX",
"49": "UT",
"50": "VT",
"51": "VA",
"53": "WA",
"54": "WV",
"55": "WI",
"56": "WY",
"60": "AS",
"66": "GU",
"69": "MP",
"72": "PR",
"78": "VI",
}
def normalize_fips(value: object, width: int) -> str:
"""Return a zero-padded FIPS code with the requested width."""
text = str(value).strip()
digits = "".join(ch for ch in text if ch.isdigit())
if not digits:
return ""
return digits.zfill(width)[-width:]
def load_counties(counties_geojson: Path) -> gpd.GeoDataFrame:
"""Load county polygons and normalize fields used downstream."""
if not counties_geojson.exists():
try:
print(
f"County GeoJSON not found at {counties_geojson}. "
f"Attempting download from {DEFAULT_COUNTIES_GEOJSON_URL}..."
)
gdf = gpd.read_file(DEFAULT_COUNTIES_GEOJSON_URL)
counties_geojson.parent.mkdir(parents=True, exist_ok=True)
# Cache the downloaded file for subsequent runs.
gdf.to_file(counties_geojson, driver="GeoJSON")
print(f"Downloaded and cached county GeoJSON to {counties_geojson}")
except Exception as exc:
raise FileNotFoundError(
f"County GeoJSON not found at {counties_geojson}, and download from "
f"{DEFAULT_COUNTIES_GEOJSON_URL} failed. Download the file manually "
"and rerun with --counties-geojson pointing to it."
) from exc
gdf = gpd.read_file(counties_geojson)
if gdf.crs is None:
gdf = gdf.set_crs("EPSG:4326")
else:
gdf = gdf.to_crs("EPSG:4326")
feature_id = None
if "id" in gdf.columns:
feature_id = gdf["id"]
elif "GEOID" in gdf.columns:
feature_id = gdf["GEOID"]
elif "GEOID10" in gdf.columns:
feature_id = gdf["GEOID10"]
elif "fips" in gdf.columns:
feature_id = gdf["fips"]
else:
raise ValueError("Unable to locate county FIPS identifier column in county polygons.")
gdf["county_fips"] = feature_id.map(lambda value: normalize_fips(value, 5))
gdf = gdf[gdf["county_fips"] != ""].copy()
if "NAME" in gdf.columns:
gdf["county_name"] = gdf["NAME"].fillna("").astype(str).str.strip()
elif "name" in gdf.columns:
gdf["county_name"] = gdf["name"].fillna("").astype(str).str.strip()
else:
gdf["county_name"] = gdf["county_fips"].map(lambda value: f"County {value}")
gdf["state_fips"] = gdf["county_fips"].str.slice(0, 2)
gdf["state"] = gdf["state_fips"].map(lambda code: STATE_FIPS_TO_ABBR.get(code, f"S{code}"))
gdf = gdf.sort_values("county_fips").reset_index(drop=True)
return gdf
+128
View File
@@ -0,0 +1,128 @@
"""Shared county zonal statistics for raster-based metrics.
Area weighting estimates how much of each raster cell lies inside a county by
rasterizing the county on a finer grid of sub-cells, then scales each cell by
the cosine of its latitude so cells count by their true surface area. Counties
that cross the 180th meridian are split so each side is read from its own small
raster window.
"""
from __future__ import annotations
import math
from typing import Dict, List
import numpy as np
import shapely
from affine import Affine
from rasterio.features import rasterize
from rasterio.windows import Window, from_bounds
from shapely.affinity import translate
from shapely.geometry import box, mapping
from shapely.geometry.base import BaseGeometry
from shapely.ops import unary_union
DEFAULT_SUBCELLS = 16
# Largest sub-cell grid rasterized for one county piece; bigger pieces use a coarser grid.
SUBCELL_BUDGET = 80_000_000
def split_at_antimeridian(geometry: BaseGeometry) -> List[BaseGeometry]:
"""Split a lon/lat geometry into pieces that each stay on one side of 180 degrees.
A county such as Aleutians West, AK has islands at both +179 and -179
degrees longitude. Its bounding box then spans nearly the whole globe, so
each side is returned as its own piece.
"""
minx, _, maxx, _ = geometry.bounds
if maxx - minx <= 180.0:
return [geometry]
positive: List[BaseGeometry] = []
negative: List[BaseGeometry] = []
for part in getattr(geometry, "geoms", [geometry]):
part_minx, _, part_maxx, _ = part.bounds
if part_maxx - part_minx > 180.0:
# One outline crosses the line: unwrap to 0..360, cut at 180, rewrap.
unwrapped = shapely.transform(
part,
lambda xy: np.column_stack((np.where(xy[:, 0] < 0, xy[:, 0] + 360.0, xy[:, 0]), xy[:, 1])),
)
positive.append(unwrapped.intersection(box(0.0, -90.0, 180.0, 90.0)))
negative.append(translate(unwrapped.intersection(box(180.0, -90.0, 360.0, 90.0)), xoff=-360.0))
elif part_minx >= 0:
positive.append(part)
else:
negative.append(part)
pieces = [unary_union(group) for group in (positive, negative) if group]
return [piece for piece in pieces if not piece.is_empty]
def geometry_window(source, geometry: BaseGeometry) -> Window:
"""Return the raster window covering a geometry, padded by one cell on each side."""
window = from_bounds(*geometry.bounds, transform=source.transform)
col_start = math.floor(window.col_off) - 1
row_start = math.floor(window.row_off) - 1
col_stop = math.ceil(window.col_off + window.width) + 1
row_stop = math.ceil(window.row_off + window.height) + 1
padded = Window(col_start, row_start, col_stop - col_start, row_stop - row_start)
return padded.intersection(Window(0, 0, source.width, source.height))
def subcells_for(shape: tuple[int, int], requested: int) -> int:
"""Return the finest sub-cell count, up to the request, that fits the budget."""
rows, cols = shape
subcells = requested
while subcells > 1 and rows * cols * subcells * subcells > SUBCELL_BUDGET:
subcells //= 2
return subcells
def cell_coverage_fractions(
shape: tuple[int, int], transform: Affine, geometry: BaseGeometry, subcells: int
) -> np.ndarray:
"""Estimate the fraction of each raster cell covered by a geometry."""
rows, cols = shape
fine = rasterize(
[mapping(geometry)],
out_shape=(rows * subcells, cols * subcells),
transform=transform * Affine.scale(1.0 / subcells),
fill=0,
default_value=1,
dtype="uint8",
)
return fine.reshape(rows, subcells, cols, subcells).mean(axis=(1, 3))
def area_weighted_class_weights(
source,
geometry: BaseGeometry,
*,
geographic: bool = True,
subcells: int = DEFAULT_SUBCELLS,
) -> Dict[int, float]:
"""Return the area inside a geometry covered by each value of a categorical raster.
Weights are relative surface areas: the fraction of each cell inside the
geometry, times cos(latitude) for geographic rasters. Cells equal to 0 or
the raster's nodata value are excluded.
"""
weights: Dict[int, float] = {}
pieces = split_at_antimeridian(geometry) if geographic else [geometry]
for piece in pieces:
window = geometry_window(source, piece)
values = source.read(1, window=window, masked=True).filled(0)
if source.nodata is not None:
values = np.where(values == source.nodata, 0, values)
transform = source.window_transform(window)
cell_weights = cell_coverage_fractions(values.shape, transform, piece, subcells_for(values.shape, subcells))
if geographic:
row_lat = transform.f + (np.arange(values.shape[0]) + 0.5) * transform.e
cell_weights = cell_weights * np.cos(np.radians(row_lat))[:, None]
counted = (values != 0) & (cell_weights > 0)
for value in np.unique(values[counted]):
code = int(value)
weights[code] = weights.get(code, 0.0) + float(cell_weights[counted & (values == value)].sum())
return weights
+64
View File
@@ -0,0 +1,64 @@
"""Koppen-Geiger raster codes and the Beck et al. legend loader."""
from __future__ import annotations
import re
from pathlib import Path
from typing import Dict
# Beck et al legend key is expected as text file, but this default handles common codes.
DEFAULT_KOPPEN_CODE_MAP = {
1: "Af",
2: "Am",
3: "Aw",
4: "BWh",
5: "BWk",
6: "BSh",
7: "BSk",
8: "Csa",
9: "Csb",
10: "Csc",
11: "Cwa",
12: "Cwb",
13: "Cwc",
14: "Cfa",
15: "Cfb",
16: "Cfc",
17: "Dsa",
18: "Dsb",
19: "Dsc",
20: "Dsd",
21: "Dwa",
22: "Dwb",
23: "Dwc",
24: "Dwd",
25: "Dfa",
26: "Dfb",
27: "Dfc",
28: "Dfd",
29: "ET",
30: "EF",
}
def load_koppen_legend(legend_path: Path | None) -> Dict[int, str]:
"""Load Koppen raster codes, using defaults when no legend exists."""
if legend_path is None:
return DEFAULT_KOPPEN_CODE_MAP
mapping: Dict[int, str] = {}
for line in legend_path.read_text(encoding="utf-8").splitlines():
text = line.strip()
if not text or text.startswith("#"):
continue
# Handles patterns like:
# "1: Af ..." or "1 = Af" or "1 Af"
match = re.match(r"^(\d+)\s*[:=]?\s*([A-Za-z]{2,3})\b", text)
if not match:
continue
key = int(match.group(1))
value = match.group(2)
mapping[key] = value
return mapping if mapping else DEFAULT_KOPPEN_CODE_MAP
+33 -18
View File
@@ -12,13 +12,24 @@ The browser blocks `fetch("data/climate-data.csv")` when `index.html` is opened
Then open [http://localhost:8000/](http://localhost:8000/). This keeps the app CSV-only while allowing the map and filters to load normally. Then open [http://localhost:8000/](http://localhost:8000/). This keeps the app CSV-only while allowing the map and filters to load normally.
## Source 1: Koppen-Geiger classes (`koppenZone`) ## Source 1: Koppen-Geiger classes (`koppenZone`, `koppenPrimaryClass`, `koppenSecondaryClass`)
- Dataset: Beck et al. updated 1-km Koppen-Geiger climate classes (historical + future windows) - Dataset: Beck et al. updated 1-km Koppen-Geiger climate classes (historical + future windows)
- Landing page: [https://www.gloh2o.org/koppen/](https://www.gloh2o.org/koppen/) - Landing page: [https://www.gloh2o.org/koppen/](https://www.gloh2o.org/koppen/)
- Primary paper for updated release: [https://www.nature.com/articles/s41597-023-02549-6](https://www.nature.com/articles/s41597-023-02549-6) - Primary paper for updated release: [https://www.nature.com/articles/s41597-023-02549-6](https://www.nature.com/articles/s41597-023-02549-6)
- Coverage: 1901-2099 (use historical 1991-2020 layer for this project to align with NOAA baselines) - Coverage: 1901-2099 (use historical 1991-2020 layer for this project to align with NOAA baselines)
- License shown on dataset page: CC BY 4.0 - License shown on dataset page: CC BY 4.0
- Local files: `data/koppen_geiger_tif/1991_2020/koppen_geiger_0p00833333.tif` and `data/koppen_geiger_tif/legend.txt`
Build the county metric and apply it to the app CSV:
```powershell
.venv\Scripts\python.exe scripts\build_county_koppen_metric.py
.venv\Scripts\python.exe scripts\apply_koppen_metric_to_climate_data.py --dry-run
.venv\Scripts\python.exe scripts\apply_koppen_metric_to_climate_data.py
```
The builder writes `data/metrics/koppen.csv` with each county's class, top and runner-up classes, and their area-weighted shares. The apply step writes `koppenZone`, `koppenPrimaryClass`, and `koppenSecondaryClass`; `--dry-run` reports the changes without writing. Run it after `build_county_climate_data.py`, which still writes an older largest-share `koppenZone`. The classification rule is documented in [`docs/filter-calculations.md`](../docs/filter-calculations.md) §1.
## Source 2: NOAA 1991-2020 gridded normals (`avgTempF`, `annualPrecipIn`, `seasonalityIndex`, previous `extremeDays`) ## Source 2: NOAA 1991-2020 gridded normals (`avgTempF`, `annualPrecipIn`, `seasonalityIndex`, previous `extremeDays`)
@@ -45,14 +56,14 @@ If you have monthly nClimGrid history files (for example `nclimgrid_tavg.nc`) ra
- Original Plotly county geometry source: [https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json](https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json) - Original Plotly county geometry source: [https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json](https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json)
- Official county geometry reference (Census TIGER/Line): [https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html](https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html) - Official county geometry reference (Census TIGER/Line): [https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html](https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html)
## Source 4: Solar resource (`avgSolarGhiKwhM2Day`) ## Source 4: Solar resource (`meanDailyGlobalHorizontalRadiationKwhM2Day`)
- Recommended dataset: NREL National Solar Radiation Database (NSRDB) - Recommended dataset: NREL National Solar Radiation Database (NSRDB)
- Data/API page: [https://developer.nrel.gov/docs/solar/nsrdb/](https://developer.nrel.gov/docs/solar/nsrdb/) - Data/API page: [https://developer.nrel.gov/docs/solar/nsrdb/](https://developer.nrel.gov/docs/solar/nsrdb/)
- Maps/geospatial data page: [https://www.nrel.gov/gis/solar-resource-maps](https://www.nrel.gov/gis/solar-resource-maps) - Maps/geospatial data page: [https://www.nrel.gov/gis/solar-resource-maps](https://www.nrel.gov/gis/solar-resource-maps)
- Fallback point API option: NASA POWER `ALLSKY_SFC_SW_DWN` (surface shortwave downwelling radiation), [https://power.larc.nasa.gov/docs/tutorials/service-data-request/api/](https://power.larc.nasa.gov/docs/tutorials/service-data-request/api/) - Fallback point API option: NASA POWER `ALLSKY_SFC_SW_DWN` (surface shortwave downwelling radiation), [https://power.larc.nasa.gov/docs/tutorials/service-data-request/api/](https://power.larc.nasa.gov/docs/tutorials/service-data-request/api/)
The app metric is designed for annual average daily global horizontal irradiance (GHI), in `kWh/m2/day`. The app metric is designed for mean daily global horizontal radiation (GHI), in `kWh/m2/day`.
For county means, use a gridded annual GHI raster and pass it to the generator with `--solar-ghi-raster`. For county means, use a gridded annual GHI raster and pass it to the generator with `--solar-ghi-raster`.
### First test: NSRDB county representative points ### First test: NSRDB county representative points
@@ -80,12 +91,12 @@ Outputs:
- `data/nrel/county_representative_points.csv`: county FIPS, name, state, latitude, and longitude. - `data/nrel/county_representative_points.csv`: county FIPS, name, state, latitude, and longitude.
- `data/nrel/representative_point_csv/`: cached raw NSRDB CSV responses by county FIPS. - `data/nrel/representative_point_csv/`: cached raw NSRDB CSV responses by county FIPS.
- `data/nrel/county_representative_point_ghi_summary.csv`: summarized `avgSolarGhiKwhM2Day` values. - `data/nrel/county_representative_point_ghi_summary.csv`: summarized `meanDailyGlobalHorizontalRadiationKwhM2Day` values.
The calculation is: The calculation is:
```text ```text
avgSolarGhiKwhM2Day = sum(hourly GHI) / 1000 / 365 meanDailyGlobalHorizontalRadiationKwhM2Day = sum(hourly GHI) / 1000 / 365
``` ```
If the raw cache is complete but the summary CSV only contains the last fetched batch, rebuild the summary from cached files without calling the API: If the raw cache is complete but the summary CSV only contains the last fetched batch, rebuild the summary from cached files without calling the API:
@@ -96,9 +107,9 @@ If the raw cache is complete but the summary CSV only contains the last fetched
This point-based workflow is easier to validate and resume than full polygon downloads, but it is an approximation of county sunlight rather than an area-weighted county mean. This point-based workflow is easier to validate and resume than full polygon downloads, but it is an approximation of county sunlight rather than an area-weighted county mean.
### First cloud-cover pass: NSRDB representative points ### First clear-sky GHI reduction pass: NSRDB representative points
For a fast cloud-cover input layer, fetch representative-point NSRDB CSVs with observed GHI, Clearsky GHI, and Cloud Type: For a fast clear-sky GHI reduction input layer, fetch representative-point NSRDB CSVs with observed GHI, Clearsky GHI, and Cloud Type:
```powershell ```powershell
.venv\Scripts\python.exe scripts\fetch_nsrdb_representative_point_cloud_metrics.py --limit 10 .venv\Scripts\python.exe scripts\fetch_nsrdb_representative_point_cloud_metrics.py --limit 10
@@ -113,18 +124,18 @@ Once the first batch looks right, fetch every county:
Outputs: Outputs:
- `data/nrel/representative_point_cloud_csv/`: cached raw NSRDB CSV responses with `ghi,clearsky_ghi,cloud_type`. - `data/nrel/representative_point_cloud_csv/`: cached raw NSRDB CSV responses with `ghi,clearsky_ghi,cloud_type`.
- `data/nrel/county_representative_point_cloud_summary.csv`: summarized cloud-cover proxy fields. - `data/nrel/county_representative_point_cloud_summary.csv`: summarized clear-sky GHI reduction fields.
- `data/nrel/county_representative_point_cloud_error_log.csv`: failed county requests with redacted API context. - `data/nrel/county_representative_point_cloud_error_log.csv`: failed county requests with redacted API context.
The primary calculation uses daylight rows where Clearsky GHI is at least 50 W/m2: The primary calculation uses daylight rows where Clearsky GHI is at least 50 W/m2:
```text ```text
cloudinessIndexPct = 1 - mean(clamped(GHI / Clearsky GHI, 0, 1)) clearSkyGhiReductionIndex = 1 - mean(clamped(GHI / Clearsky GHI, 0, 1))
``` ```
The same summary also stores daylight row counts, observed-to-clear-sky ratio, and broad Cloud Type frequency buckets. This is faster than the polygon archive workflow, but it remains a representative-point county approximation. The same summary also stores daylight row counts, observed-to-clear-sky ratio, and broad Cloud Type frequency buckets. This is faster than the polygon archive workflow, but it remains a representative-point county approximation.
Apply the representative-point cloudiness metric to the browser app CSV: Apply the representative-point clear-sky GHI reduction metric to the browser app CSV:
```powershell ```powershell
.venv\Scripts\python.exe scripts\apply_nsrdb_cloud_metric_to_climate_data.py .venv\Scripts\python.exe scripts\apply_nsrdb_cloud_metric_to_climate_data.py
@@ -134,7 +145,7 @@ The current cloud summary covers the 3,143 county representative points and leav
### County-average target: NSRDB polygon cloud archive requests ### County-average target: NSRDB polygon cloud archive requests
For a less noisy county-level cloudiness layer, submit county polygons to the NSRDB archive workflow with `ghi,clearsky_ghi,cloud_type`. This uses the same polygon tiling and pacing logic as the GHI archive workflow, but writes separate cloud manifests, response JSON files, archives, and summaries. For a less noisy county-level clear-sky GHI reduction layer, submit county polygons to the NSRDB archive workflow with `ghi,clearsky_ghi,cloud_type`. This uses the same polygon tiling and pacing logic as the GHI archive workflow, but writes separate cloud manifests, response JSON files, archives, and summaries.
Start with a dry run: Start with a dry run:
@@ -166,7 +177,7 @@ Download completed archives:
.venv\Scripts\python.exe scripts\download_nsrdb_county_polygon_cloud_archives.py .venv\Scripts\python.exe scripts\download_nsrdb_county_polygon_cloud_archives.py
``` ```
Summarize downloaded archives into area-weighted county cloudiness: Summarize downloaded archives into area-weighted county clear-sky GHI reduction:
```powershell ```powershell
.venv\Scripts\python.exe scripts\summarize_nsrdb_county_polygon_cloud_archives.py --reuse-existing-output .venv\Scripts\python.exe scripts\summarize_nsrdb_county_polygon_cloud_archives.py --reuse-existing-output
@@ -178,7 +189,7 @@ Outputs:
- `data/nrel/county_polygon_cloud_error_log.csv`: cloud archive request errors. - `data/nrel/county_polygon_cloud_error_log.csv`: cloud archive request errors.
- `data/nrel/polygon_cloud_request_responses/`: raw NSRDB cloud acknowledgement JSON files. - `data/nrel/polygon_cloud_request_responses/`: raw NSRDB cloud acknowledgement JSON files.
- `data/nrel/polygon_cloud_archives/`: downloaded cloud ZIP archives. - `data/nrel/polygon_cloud_archives/`: downloaded cloud ZIP archives.
- `data/nrel/county_polygon_cloud_summary.csv`: area-weighted county cloudiness summaries. - `data/nrel/county_polygon_cloud_summary.csv`: area-weighted county clear-sky GHI reduction summaries.
Then update the browser app CSV. Polygon area-weighted values are used first when `data/nrel/county_polygon_cloud_summary.csv` exists; representative-point values remain the fallback: Then update the browser app CSV. Polygon area-weighted values are used first when `data/nrel/county_polygon_cloud_summary.csv` exists; representative-point values remain the fallback:
@@ -228,7 +239,7 @@ Outputs:
- `data/nrel/county_polygon_ghi_error_log.csv`: counties that need retrying or tiling. - `data/nrel/county_polygon_ghi_error_log.csv`: counties that need retrying or tiling.
- `data/nrel/polygon_request_responses/`: raw API acknowledgement JSON files. - `data/nrel/polygon_request_responses/`: raw API acknowledgement JSON files.
- `data/nrel/polygon_archives/`: downloaded county or tiled county ZIP archives. - `data/nrel/polygon_archives/`: downloaded county or tiled county ZIP archives.
- `data/nrel/county_polygon_ghi_summary.csv`: summarized polygon archive GHI, including `avgSolarGhiKwhM2Day`. - `data/nrel/county_polygon_ghi_summary.csv`: summarized polygon archive GHI, including `meanDailyGlobalHorizontalRadiationKwhM2Day`.
Download completed GHI archives: Download completed GHI archives:
@@ -255,15 +266,18 @@ Then update the app CSV. Polygon archive GHI is used first; representative-point
## Metric definitions in generated output ## Metric definitions in generated output
- `koppenZone`: majority class within county polygon from Koppen raster. - `koppenZone`: the county's predominant Koppen-Geiger class, meaning the class covering at least 50% of the county's land area and leading the runner-up by at least 5 percentage points; otherwise `Mixed`. Shares are area-weighted, with ocean and no-data cells excluded.
- `koppenPrimaryClass` / `koppenSecondaryClass`: for Mixed counties only, the top and runner-up classes, drawn as stripes on the map; blank for predominant counties.
- `avgTempF`: mean of 12 monthly county mean temperatures, converted C -> F. - `avgTempF`: mean of 12 monthly county mean temperatures, converted C -> F.
- `annualPrecipIn`: sum of 12 monthly county mean precipitation totals, converted mm -> inches. - `annualPrecipIn`: sum of 12 monthly county mean precipitation totals, converted mm -> inches.
- `seasonalityIndex`: coefficient of variation of monthly precipitation totals, scaled to 0-100. - `seasonalityIndex`: coefficient of variation of monthly precipitation totals, scaled to 0-100.
- `wettestPrecipMonth`: month with the highest 1991-2020 county mean precipitation total. - `wettestPrecipMonth`: month with the highest 1991-2020 county mean precipitation total.
- `driestPrecipMonth`: month with the lowest 1991-2020 county mean precipitation total. - `driestPrecipMonth`: month with the lowest 1991-2020 county mean precipitation total.
- previous `extremeDays` / `oldExtremeDays`: count of daily-normal or monthly-proxy days where county mean `tmax >= 95F` or `tmin <= 32F` (thresholds configurable in script). This is preserved only as an audit column after the NOAA nClimGrid-Daily metrics are applied. - previous `extremeDays` / `oldExtremeDays`: count of daily-normal or monthly-proxy days where county mean `tmax >= 95F` or `tmin <= 32F` (thresholds configurable in script). This is preserved only as an audit column after the NOAA nClimGrid-Daily metrics are applied.
- `avgSolarGhiKwhM2Day`: county mean annual average daily GHI, in `kWh/m2/day`. The current app CSV uses NSRDB polygon archive area-weighted values where available, with representative-point values kept as fallback. - `humidHeatDays`: average annual count of days where estimated Heat Index is at least 90 F, the lower bound of the NWS Extreme Caution category. The NWS Heat Index algorithm is applied to NOAA nClimGrid-Daily county `tmax` and gridMET daily minimum relative humidity (`rmin`) for 1991-2020. Because the inputs are paired daily extrema rather than coincident hourly observations, this is an estimated daily-peak proxy.
- `cloudinessIndexPct`: representative-point NSRDB daylight cloudiness proxy derived from observed GHI divided by Clearsky GHI. Higher values mean observed irradiance is lower relative to modeled clear-sky irradiance. Despite the legacy field name, values are stored on a 0-1 scale. - `avgSummerSpecificHumidityGKg`: county-cell-weighted mean gridMET specific humidity for June-August 1991-2020, converted from kg/kg to g/kg.
- `meanDailyGlobalHorizontalRadiationKwhM2Day`: county mean daily global horizontal radiation (GHI), in `kWh/m2/day`. The current app CSV uses NSRDB polygon archive area-weighted values where available, with representative-point values kept as fallback.
- `clearSkyGhiReductionIndex`: NSRDB daylight clear-sky GHI reduction index derived from observed GHI divided by modeled clear-sky GHI. Higher values mean observed irradiance is lower relative to clear-sky conditions. Values are stored on a 0-1 scale.
## Run the generator ## Run the generator
@@ -293,11 +307,12 @@ python scripts/build_county_climate_data.py `
Notes: Notes:
- This writes only the base columns and an older largest-share `koppenZone`. Do not run it over the live `data/climate-data.csv`: it would drop the columns added by later stages. After a full rebuild, run the Köppen build and apply steps (Source 1) and the enrichment stages.
- This computes all counties in your geometry file, not just the sample records. - This computes all counties in your geometry file, not just the sample records.
- For counties outside CONUS coverage in NOAA gridded files, fallback values are applied by the script when no valid grid values intersect. - For counties outside CONUS coverage in NOAA gridded files, fallback values are applied by the script when no valid grid values intersect.
- For physically-based daily `extremeDays`, provide true daily grids and set `--extreme-days-mode require-daily`. - For physically-based daily `extremeDays`, provide true daily grids and set `--extreme-days-mode require-daily`.
- If `--counties-geojson` does not exist locally, the script will try to download the county GeoJSON automatically from the Plotly URL above and cache it at that path. - If `--counties-geojson` does not exist locally, the script will try to download the county GeoJSON automatically from the Plotly URL above and cache it at that path.
- `--solar-ghi-csv` is optional and can load county-keyed solar summaries into `avgSolarGhiKwhM2Day`. - `--solar-ghi-csv` is optional and can load county-keyed solar summaries into `meanDailyGlobalHorizontalRadiationKwhM2Day`.
- `--solar-ghi-raster` is optional and takes precedence over `--solar-ghi-csv`. If you have a gridded annual GHI raster, add `--solar-ghi-raster path/to/annual_ghi_kwh_m2_day.tif` for a true county-area raster mean. - `--solar-ghi-raster` is optional and takes precedence over `--solar-ghi-csv`. If you have a gridded annual GHI raster, add `--solar-ghi-raster path/to/annual_ghi_kwh_m2_day.tif` for a true county-area raster mean.
## Source 5: NOAA nClimGrid-Daily county area averages (`locallyExtremeDays`) ## Source 5: NOAA nClimGrid-Daily county area averages (`locallyExtremeDays`)
-1
View File
@@ -20,7 +20,6 @@ import urllib.error
import urllib.request import urllib.request
from pathlib import Path from pathlib import Path
DEFAULT_BASE_URL = "https://www.northwestknowledge.net/metdata/data" DEFAULT_BASE_URL = "https://www.northwestknowledge.net/metdata/data"
DEFAULT_VARIABLES = ("sph", "rmax", "rmin") DEFAULT_VARIABLES = ("sph", "rmax", "rmin")
@@ -20,7 +20,6 @@ from enum import Enum
from pathlib import Path from pathlib import Path
from typing import Any from typing import Any
DEFAULT_TIMEOUT = 300 DEFAULT_TIMEOUT = 300
DEFAULT_CHUNK_SIZE = 1024 * 1024 DEFAULT_CHUNK_SIZE = 1024 * 1024
@@ -12,7 +12,6 @@ import sys
import download_nsrdb_county_polygon_archives as polygon_download import download_nsrdb_county_polygon_archives as polygon_download
DEFAULT_ARGS = [ DEFAULT_ARGS = [
"--response-dir", "--response-dir",
"data/nrel/polygon_cloud_request_responses", "data/nrel/polygon_cloud_request_responses",
@@ -13,7 +13,6 @@ import sys
import download_nsrdb_county_polygon_archives as polygon_download import download_nsrdb_county_polygon_archives as polygon_download
DEFAULT_ARGS = [ DEFAULT_ARGS = [
"--response-dir", "--response-dir",
"data/nrel/polygon_request_responses", "data/nrel/polygon_request_responses",
@@ -1,6 +1,6 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
""" """
Fetch NSRDB representative-point inputs for a county cloud-cover metric. Fetch NSRDB representative-point inputs for a county clear-sky GHI reduction metric.
This uses the direct single-point NSRDB CSV endpoint rather than polygon archive This uses the direct single-point NSRDB CSV endpoint rather than polygon archive
requests. It is faster for first-pass cloud metrics because it downloads one CSV requests. It is faster for first-pass cloud metrics because it downloads one CSV
@@ -12,10 +12,10 @@ ghi,clearsky_ghi,cloud_type
The summary metric is based on daylight rows with valid GHI and Clearsky GHI: The summary metric is based on daylight rows with valid GHI and Clearsky GHI:
cloudinessIndexPct = 1 - mean(clamped(GHI / Clearsky GHI)) clearSkyGhiReductionIndex = 1 - mean(clamped(GHI / Clearsky GHI))
where the ratio is clamped to [0, 1] so occasional above-clear-sky modeled GHI where the ratio is clamped to [0, 1] so occasional above-clear-sky modeled GHI
does not create negative cloudiness. does not create negative reduction values.
""" """
from __future__ import annotations from __future__ import annotations
@@ -33,7 +33,6 @@ import urllib.parse
import urllib.request import urllib.request
from pathlib import Path from pathlib import Path
DEFAULT_POINTS_CSV = Path("data/nrel/county_representative_points.csv") DEFAULT_POINTS_CSV = Path("data/nrel/county_representative_points.csv")
DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_point_cloud_summary.csv") DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_point_cloud_summary.csv")
DEFAULT_ERROR_CSV = Path("data/nrel/county_representative_point_cloud_error_log.csv") DEFAULT_ERROR_CSV = Path("data/nrel/county_representative_point_cloud_error_log.csv")
@@ -73,7 +72,7 @@ CLOUD_SUMMARY_FIELDS = [
"state_abbr", "state_abbr",
"lat", "lat",
"lon", "lon",
"cloudinessIndexPct", "clearSkyGhiReductionIndex",
"avgObservedToClearskyRatio", "avgObservedToClearskyRatio",
"daylightRows", "daylightRows",
"allRows", "allRows",
@@ -442,7 +441,7 @@ def summarize_cloud_metrics(
raise ValueError("No daylight rows with valid GHI and Clearsky GHI were found.") raise ValueError("No daylight rows with valid GHI and Clearsky GHI were found.")
avg_ratio = ratio_sum / daylight_rows avg_ratio = ratio_sum / daylight_rows
cloudiness_pct = 1 - avg_ratio clear_sky_ghi_reduction = 1 - avg_ratio
clear_or_probably_clear = daylight_cloud_type_counts.get("0", 0) + daylight_cloud_type_counts.get("1", 0) clear_or_probably_clear = daylight_cloud_type_counts.get("0", 0) + daylight_cloud_type_counts.get("1", 0)
cloudy_or_obscured = sum( cloudy_or_obscured = sum(
daylight_cloud_type_counts.get(code, 0) daylight_cloud_type_counts.get(code, 0)
@@ -467,7 +466,7 @@ def summarize_cloud_metrics(
"state_abbr": point["state_abbr"], "state_abbr": point["state_abbr"],
"lat": point["lat"], "lat": point["lat"],
"lon": point["lon"], "lon": point["lon"],
"cloudinessIndexPct": fmt(cloudiness_pct, 4), "clearSkyGhiReductionIndex": fmt(clear_sky_ghi_reduction, 4),
"avgObservedToClearskyRatio": fmt(avg_ratio, 4), "avgObservedToClearskyRatio": fmt(avg_ratio, 4),
"daylightRows": str(daylight_rows), "daylightRows": str(daylight_rows),
"allRows": str(all_rows), "allRows": str(all_rows),
@@ -675,7 +674,7 @@ def log_county_result(
if summary is not None: if summary is not None:
print( print(
" " " "
f"cloudiness={summary['cloudinessIndexPct']}, " f"clear_sky_ghi_reduction={summary['clearSkyGhiReductionIndex']}, "
f"clear_ratio={summary['avgObservedToClearskyRatio']}, " f"clear_ratio={summary['avgObservedToClearskyRatio']}, "
f"daylight_rows={summary['daylightRows']}" f"daylight_rows={summary['daylightRows']}"
) )
@@ -783,7 +782,7 @@ def parse_args() -> argparse.Namespace:
"--min-clearsky-ghi", "--min-clearsky-ghi",
type=float, type=float,
default=DEFAULT_MIN_CLEARSKY_GHI, default=DEFAULT_MIN_CLEARSKY_GHI,
help="Minimum Clearsky GHI W/m2 for daylight cloudiness ratio rows.", help="Minimum Clearsky GHI W/m2 for daylight GHI ratio rows.",
) )
parser.add_argument( parser.add_argument(
"--no-polar-fallback", "--no-polar-fallback",
@@ -20,7 +20,6 @@ import urllib.parse
import urllib.request import urllib.request
from pathlib import Path from pathlib import Path
DEFAULT_POINTS_CSV = Path("data/nrel/county_representative_points.csv") DEFAULT_POINTS_CSV = Path("data/nrel/county_representative_points.csv")
DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_point_ghi_summary.csv") DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_point_ghi_summary.csv")
DEFAULT_ERROR_CSV = Path("data/nrel/county_representative_point_ghi_error_log.csv") DEFAULT_ERROR_CSV = Path("data/nrel/county_representative_point_ghi_error_log.csv")
@@ -247,7 +246,7 @@ def summarize_ghi(point: dict[str, str], csv_text: str, source_file: Path, sourc
"state_abbr": point["state_abbr"], "state_abbr": point["state_abbr"],
"lat": point["lat"], "lat": point["lat"],
"lon": point["lon"], "lon": point["lon"],
"avgSolarGhiKwhM2Day": f"{avg_daily_ghi:.3f}", "meanDailyGlobalHorizontalRadiationKwhM2Day": f"{avg_daily_ghi:.3f}",
"ghi_rows": str(len(ghi_values)), "ghi_rows": str(len(ghi_values)),
"ghi_min": f"{min(ghi_values):.1f}", "ghi_min": f"{min(ghi_values):.1f}",
"ghi_max": f"{max(ghi_values):.1f}", "ghi_max": f"{max(ghi_values):.1f}",
@@ -269,7 +268,7 @@ def write_summary_csv(rows: list[dict[str, str]], output_csv: Path) -> None:
"state_abbr", "state_abbr",
"lat", "lat",
"lon", "lon",
"avgSolarGhiKwhM2Day", "meanDailyGlobalHorizontalRadiationKwhM2Day",
"ghi_rows", "ghi_rows",
"ghi_min", "ghi_min",
"ghi_max", "ghi_max",
@@ -4,7 +4,7 @@ Rebuild county representative-point GHI summaries from cached NSRDB raw CSV resp
This does not call the NSRDB API. It reads data/nrel/representative_point_csv/*.csv, computes: This does not call the NSRDB API. It reads data/nrel/representative_point_csv/*.csv, computes:
avgSolarGhiKwhM2Day = sum(hourly GHI) / 1000 / 365 meanDailyGlobalHorizontalRadiationKwhM2Day = sum(hourly GHI) / 1000 / 365
and writes a county-keyed CSV compatible with build_county_climate_data.py. and writes a county-keyed CSV compatible with build_county_climate_data.py.
""" """
@@ -15,7 +15,6 @@ import argparse
import csv import csv
from pathlib import Path from pathlib import Path
DEFAULT_POINTS_CSV = Path("data/nrel/county_representative_points.csv") DEFAULT_POINTS_CSV = Path("data/nrel/county_representative_points.csv")
DEFAULT_REPRESENTATIVE_POINT_CSV_DIR = Path("data/nrel/representative_point_csv") DEFAULT_REPRESENTATIVE_POINT_CSV_DIR = Path("data/nrel/representative_point_csv")
DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_point_ghi_summary.csv") DEFAULT_OUTPUT_CSV = Path("data/nrel/county_representative_point_ghi_summary.csv")
@@ -82,7 +81,7 @@ def summarize_county(point: dict[str, str], raw_csv: Path) -> dict[str, str]:
"state_abbr": point["state_abbr"], "state_abbr": point["state_abbr"],
"lat": point["lat"], "lat": point["lat"],
"lon": point["lon"], "lon": point["lon"],
"avgSolarGhiKwhM2Day": f"{avg_daily_ghi:.3f}", "meanDailyGlobalHorizontalRadiationKwhM2Day": f"{avg_daily_ghi:.3f}",
"ghi_rows": str(len(ghi_values)), "ghi_rows": str(len(ghi_values)),
"ghi_min": f"{min(ghi_values):.1f}", "ghi_min": f"{min(ghi_values):.1f}",
"ghi_max": f"{max(ghi_values):.1f}", "ghi_max": f"{max(ghi_values):.1f}",
@@ -104,7 +103,7 @@ def write_summary_csv(rows: list[dict[str, str]], output_csv: Path) -> None:
"state_abbr", "state_abbr",
"lat", "lat",
"lon", "lon",
"avgSolarGhiKwhM2Day", "meanDailyGlobalHorizontalRadiationKwhM2Day",
"ghi_rows", "ghi_rows",
"ghi_min", "ghi_min",
"ghi_max", "ghi_max",
@@ -29,7 +29,6 @@ import geopandas as gpd
from shapely.geometry import box from shapely.geometry import box
from shapely.wkt import dumps as dump_wkt from shapely.wkt import dumps as dump_wkt
DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json") DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json")
DEFAULT_ENDPOINT = "https://developer.nlr.gov/api/nsrdb/v2/solar/nsrdb-GOES-tmy-v4-0-0-download.json" DEFAULT_ENDPOINT = "https://developer.nlr.gov/api/nsrdb/v2/solar/nsrdb-GOES-tmy-v4-0-0-download.json"
DEFAULT_POLAR_ENDPOINT = "https://developer.nlr.gov/api/nsrdb/v2/solar/nsrdb-polar-tmy-v4-0-0-download.json" DEFAULT_POLAR_ENDPOINT = "https://developer.nlr.gov/api/nsrdb/v2/solar/nsrdb-polar-tmy-v4-0-0-download.json"
@@ -550,7 +549,7 @@ class ArchiveQueueMonitor:
) )
) )
now = time.monotonic() now = time.monotonic()
for url, status in zip(unresolved_urls, statuses): for url, status in zip(unresolved_urls, statuses, strict=True):
if status == "pending": if status == "pending":
if url not in self.pending_urls: if url not in self.pending_urls:
ambiguous_403 += 1 ambiguous_403 += 1
@@ -584,7 +583,7 @@ class ArchiveQueueMonitor:
return { return {
"pending": sum( "pending": sum(
status == "pending" and url in self.pending_urls status == "pending" and url in self.pending_urls
for url, status in zip(unresolved_urls, statuses) for url, status in zip(unresolved_urls, statuses, strict=True)
), ),
"ambiguous_403": ambiguous_403, "ambiguous_403": ambiguous_403,
"recently_downloaded": recently_downloaded, "recently_downloaded": recently_downloaded,
@@ -1402,8 +1401,6 @@ def submit_county_request(
return [], build_error_row(county, profile, None, None, wkt, error) return [], build_error_row(county, profile, None, None, wkt, error)
for profile in request_profiles(args, county): for profile in request_profiles(args, county):
site_count: int | None = None
weight: int | None = None
final_profile = profile final_profile = profile
try: try:
return [ return [
@@ -1,6 +1,6 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
""" """
Submit NSRDB polygon archive requests for county-average cloudiness inputs. Submit NSRDB polygon archive requests for county-average clear-sky GHI reduction inputs.
This is a cloud-specific wrapper around request_nsrdb_county_polygon_archives.py. This is a cloud-specific wrapper around request_nsrdb_county_polygon_archives.py.
It reuses the existing polygon tiling, site-count, Polar fallback, pacing, and It reuses the existing polygon tiling, site-count, Polar fallback, pacing, and
@@ -15,7 +15,6 @@ import sys
import request_nsrdb_county_polygon_archives as polygon_request import request_nsrdb_county_polygon_archives as polygon_request
DEFAULT_ARGS = [ DEFAULT_ARGS = [
"--requests-csv", "--requests-csv",
"data/nrel/county_polygon_cloud_request_manifest.csv", "data/nrel/county_polygon_cloud_request_manifest.csv",
@@ -13,7 +13,6 @@ import sys
import request_nsrdb_county_polygon_archives as polygon_request import request_nsrdb_county_polygon_archives as polygon_request
DEFAULT_ARGS = [ DEFAULT_ARGS = [
"--requests-csv", "--requests-csv",
"data/nrel/county_polygon_ghi_request_manifest.csv", "data/nrel/county_polygon_ghi_request_manifest.csv",
+28 -7
View File
@@ -25,7 +25,6 @@ import numpy as np
import pandas as pd import pandas as pd
import xarray as xr import xarray as xr
STATE_FIPS_TO_ABBR = { STATE_FIPS_TO_ABBR = {
"01": "AL", "01": "AL",
"04": "AZ", "04": "AZ",
@@ -80,6 +79,7 @@ STATE_FIPS_TO_ABBR = {
CONUS_STATE_FIPS = set(STATE_FIPS_TO_ABBR) CONUS_STATE_FIPS = set(STATE_FIPS_TO_ABBR)
REQUIRED_VARIABLES = ("sph", "rmax", "rmin") REQUIRED_VARIABLES = ("sph", "rmax", "rmin")
SUMMER_MONTHS = {6, 7, 8} SUMMER_MONTHS = {6, 7, 8}
DEFAULT_HEAT_INDEX_THRESHOLD_F = 90.0
NOAA_REGION_CODE_TO_FIPS_OVERRIDES = { NOAA_REGION_CODE_TO_FIPS_OVERRIDES = {
"18511": "11001", "18511": "11001",
} }
@@ -436,8 +436,12 @@ def _county_means_chunk(
def _heat_index_f(t_f: np.ndarray, rh_pct: np.ndarray) -> np.ndarray: def _heat_index_f(t_f: np.ndarray, rh_pct: np.ndarray) -> np.ndarray:
"""Return the NWS Heat Index for Fahrenheit temperature and percent RH."""
rh = np.clip(rh_pct, 0, 100) rh = np.clip(rh_pct, 0, 100)
heat_index = ( simple_heat_index = 0.5 * (t_f + 61.0 + ((t_f - 68.0) * 1.2) + (rh * 0.094))
simple_heat_index = (simple_heat_index + t_f) / 2
regression_heat_index = (
-42.379 -42.379
+ 2.04901523 * t_f + 2.04901523 * t_f
+ 10.14333127 * rh + 10.14333127 * rh
@@ -451,9 +455,17 @@ def _heat_index_f(t_f: np.ndarray, rh_pct: np.ndarray) -> np.ndarray:
low_rh_adjustment = ((13 - rh) / 4) * np.sqrt(np.maximum((17 - np.abs(t_f - 95)) / 17, 0)) low_rh_adjustment = ((13 - rh) / 4) * np.sqrt(np.maximum((17 - np.abs(t_f - 95)) / 17, 0))
high_rh_adjustment = ((rh - 85) / 10) * ((87 - t_f) / 5) high_rh_adjustment = ((rh - 85) / 10) * ((87 - t_f) / 5)
heat_index = np.where((rh < 13) & (80 <= t_f) & (t_f <= 112), heat_index - low_rh_adjustment, heat_index) regression_heat_index = np.where(
heat_index = np.where((rh > 85) & (80 <= t_f) & (t_f <= 87), heat_index + high_rh_adjustment, heat_index) (rh < 13) & (80 <= t_f) & (t_f <= 112),
return np.where(t_f >= 80, heat_index, t_f) regression_heat_index - low_rh_adjustment,
regression_heat_index,
)
regression_heat_index = np.where(
(rh > 85) & (80 <= t_f) & (t_f <= 87),
regression_heat_index + high_rh_adjustment,
regression_heat_index,
)
return np.where(simple_heat_index >= 80, regression_heat_index, simple_heat_index)
def _time_dimension(data_array: xr.DataArray, lat_dim: str, lon_dim: str) -> str: def _time_dimension(data_array: xr.DataArray, lat_dim: str, lon_dim: str) -> str:
@@ -646,7 +658,8 @@ def summarize(args: argparse.Namespace) -> None:
"gridmetSource": ( "gridmetSource": (
f"gridmet-{source_period} county-cell-weighted; " f"gridmet-{source_period} county-cell-weighted; "
"relative-humidity=(rmax+rmin)/2; " "relative-humidity=(rmax+rmin)/2; "
"humidHeatDays=noaa-nclimgrid-tmax+gridmet-rmin heat-index-gte-90f" "humidHeatDays=noaa-nclimgrid-tmax+gridmet-rmin "
f"heat-index-gte-{args.heat_index_threshold_f:g}f"
), ),
} }
) )
@@ -680,7 +693,15 @@ def main() -> int:
parser.add_argument("--years", type=_parse_years, default=None, help="Optional comma/range years, e.g. 1991-2020.") parser.add_argument("--years", type=_parse_years, default=None, help="Optional comma/range years, e.g. 1991-2020.")
parser.add_argument("--start-year", type=int, default=None) parser.add_argument("--start-year", type=int, default=None)
parser.add_argument("--end-year", type=int, default=None) parser.add_argument("--end-year", type=int, default=None)
parser.add_argument("--heat-index-threshold-f", type=float, default=90.0) parser.add_argument(
"--heat-index-threshold-f",
type=float,
default=DEFAULT_HEAT_INDEX_THRESHOLD_F,
help=(
"Daily Heat Index threshold in F. The default 90 F is the lower "
"bound of the NWS Extreme Caution category."
),
)
parser.add_argument("--chunk-days", type=int, default=31) parser.add_argument("--chunk-days", type=int, default=31)
summarize(parser.parse_args()) summarize(parser.parse_args())
return 0 return 0
@@ -4,7 +4,7 @@ Summarize downloaded NSRDB county polygon GHI archives.
Each downloaded polygon archive contains one CSV per NSRDB site that intersects Each downloaded polygon archive contains one CSV per NSRDB site that intersects
the county polygon or tile. This script computes area-weighted county-level the county polygon or tile. This script computes area-weighted county-level
average daily GHI values across those site CSVs, combining multiple tile mean daily GHI values across those site CSVs, combining multiple tile
archives back into one county summary when present. archives back into one county summary when present.
""" """
@@ -18,7 +18,6 @@ import statistics
import zipfile import zipfile
from pathlib import Path from pathlib import Path
DEFAULT_ARCHIVE_DIR = Path("data/nrel/polygon_archives") DEFAULT_ARCHIVE_DIR = Path("data/nrel/polygon_archives")
DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json") DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json")
DEFAULT_REQUESTS_CSV = Path("data/nrel/county_polygon_ghi_request_manifest.csv") DEFAULT_REQUESTS_CSV = Path("data/nrel/county_polygon_ghi_request_manifest.csv")
@@ -37,7 +36,7 @@ SUMMARY_FIELDS = [
"request_site_count", "request_site_count",
"ghi_rows_per_site_min", "ghi_rows_per_site_min",
"ghi_rows_per_site_max", "ghi_rows_per_site_max",
"avgSolarGhiKwhM2Day", "meanDailyGlobalHorizontalRadiationKwhM2Day",
"areaWeightedAvgSolarGhiKwhM2Day", "areaWeightedAvgSolarGhiKwhM2Day",
"areaWeightedSites", "areaWeightedSites",
"weightedCellAreaKm2", "weightedCellAreaKm2",
@@ -144,7 +143,10 @@ def extract_site_lon_lat(csv_text: str) -> tuple[float, float]:
if len(rows) < 2: if len(rows) < 2:
raise ValueError("Could not read NSRDB metadata rows.") raise ValueError("Could not read NSRDB metadata rows.")
metadata = {key.strip().lower(): value.strip() for key, value in zip(rows[0], rows[1])} metadata = {
key.strip().lower(): value.strip()
for key, value in zip(rows[0], rows[1], strict=True)
}
try: try:
lon = float(metadata["longitude"]) lon = float(metadata["longitude"])
lat = float(metadata["latitude"]) lat = float(metadata["latitude"])
@@ -185,7 +187,7 @@ def extract_ghi_values(csv_text: str) -> list[float]:
def site_average_daily_ghi(csv_text: str) -> tuple[float, int]: def site_average_daily_ghi(csv_text: str) -> tuple[float, int]:
"""Return average daily GHI in kWh/m2/day and number of GHI rows.""" """Return mean daily GHI in kWh/m2/day and number of GHI rows."""
ghi_values = extract_ghi_values(csv_text) ghi_values = extract_ghi_values(csv_text)
return sum(ghi_values) / 1000 / 365, len(ghi_values) return sum(ghi_values) / 1000 / 365, len(ghi_values)
@@ -353,7 +355,7 @@ def area_weighted_average(
weighted_sites = 0 weighted_sites = 0
county_area = sum(projected_polygon_area(polygon) for polygon in projected_county) county_area = sum(projected_polygon_area(polygon) for polygon in projected_county)
for site_average, (lon, lat) in zip(site_averages, site_lon_lats): for site_average, (lon, lat) in zip(site_averages, site_lon_lats, strict=True):
x, y = project_lon_lat(lon, lat, reference_lat) x, y = project_lon_lat(lon, lat, reference_lat)
min_x = x - half_cell min_x = x - half_cell
min_y = y - half_cell min_y = y - half_cell
@@ -408,7 +410,7 @@ def summarize_archives(
"polygon_sites": len(site_averages), "polygon_sites": len(site_averages),
"ghi_rows_per_site_min": min(row_counts), "ghi_rows_per_site_min": min(row_counts),
"ghi_rows_per_site_max": max(row_counts), "ghi_rows_per_site_max": max(row_counts),
"avgSolarGhiKwhM2Day": weighted["area_weighted_avg"], "meanDailyGlobalHorizontalRadiationKwhM2Day": weighted["area_weighted_avg"],
"area_weighted_avg": weighted["area_weighted_avg"], "area_weighted_avg": weighted["area_weighted_avg"],
"area_weighted_sites": weighted["area_weighted_sites"], "area_weighted_sites": weighted["area_weighted_sites"],
"weighted_cell_area_km2": weighted["weighted_cell_area_km2"], "weighted_cell_area_km2": weighted["weighted_cell_area_km2"],
@@ -450,11 +452,11 @@ def build_summary_row(
if str(row.get("site_count", "")).strip().isdigit() if str(row.get("site_count", "")).strip().isdigit()
] ]
archive_zip = ";".join(str(path) for path in archive_paths) archive_zip = ";".join(str(path) for path in archive_paths)
final_avg = float(polygon_summary["avgSolarGhiKwhM2Day"]) final_avg = float(polygon_summary["meanDailyGlobalHorizontalRadiationKwhM2Day"])
area_weighted_avg = float(polygon_summary["area_weighted_avg"]) area_weighted_avg = float(polygon_summary["area_weighted_avg"])
representative_point_avg = ( representative_point_avg = (
float(representative_point_row["avgSolarGhiKwhM2Day"]) float(representative_point_row["meanDailyGlobalHorizontalRadiationKwhM2Day"])
if representative_point_row.get("avgSolarGhiKwhM2Day") if representative_point_row.get("meanDailyGlobalHorizontalRadiationKwhM2Day")
else None else None
) )
weighted_minus_representative_point = ( weighted_minus_representative_point = (
@@ -475,7 +477,7 @@ def build_summary_row(
"request_site_count": str(sum(request_site_counts)) if request_site_counts else request_row.get("site_count", ""), "request_site_count": str(sum(request_site_counts)) if request_site_counts else request_row.get("site_count", ""),
"ghi_rows_per_site_min": str(polygon_summary["ghi_rows_per_site_min"]), "ghi_rows_per_site_min": str(polygon_summary["ghi_rows_per_site_min"]),
"ghi_rows_per_site_max": str(polygon_summary["ghi_rows_per_site_max"]), "ghi_rows_per_site_max": str(polygon_summary["ghi_rows_per_site_max"]),
"avgSolarGhiKwhM2Day": format_float(final_avg), "meanDailyGlobalHorizontalRadiationKwhM2Day": format_float(final_avg),
"areaWeightedAvgSolarGhiKwhM2Day": format_float(area_weighted_avg), "areaWeightedAvgSolarGhiKwhM2Day": format_float(area_weighted_avg),
"areaWeightedSites": str(polygon_summary["area_weighted_sites"]), "areaWeightedSites": str(polygon_summary["area_weighted_sites"]),
"weightedCellAreaKm2": format_float(float(polygon_summary["weighted_cell_area_km2"]), 1), "weightedCellAreaKm2": format_float(float(polygon_summary["weighted_cell_area_km2"]), 1),
@@ -4,15 +4,15 @@ Summarize downloaded NSRDB county polygon cloud archives.
Each downloaded polygon archive contains one CSV per NSRDB grid site that Each downloaded polygon archive contains one CSV per NSRDB grid site that
intersects the county polygon or tile. This script computes area-weighted intersects the county polygon or tile. This script computes area-weighted
county-level cloudiness metrics across those site CSVs, combining multiple tile county-level clear-sky GHI reduction metrics across those site CSVs, combining multiple tile
archives back into one county summary when present. archives back into one county summary when present.
Primary metric: Primary metric:
cloudinessIndexPct = 1 - mean(clamped(GHI / Clearsky GHI, 0, 1)) clearSkyGhiReductionIndex = 1 - mean(clamped(GHI / Clearsky GHI, 0, 1))
Despite the legacy "Pct" field name, the primary index is stored on a 0-1 scale. The primary index is stored on a 0-1 scale. Cloud Type bucket fields are stored
Cloud Type bucket fields are stored as percentages. as percentages.
""" """
from __future__ import annotations from __future__ import annotations
@@ -26,9 +26,9 @@ import zipfile
from pathlib import Path from pathlib import Path
from summarize_nsrdb_county_polygon_archives import ( from summarize_nsrdb_county_polygon_archives import (
area_weighted_average,
archive_groups_by_county, archive_groups_by_county,
archive_signature, archive_signature,
area_weighted_average,
county_fips_from_archive, county_fips_from_archive,
existing_archive_signature, existing_archive_signature,
extract_site_lon_lat, extract_site_lon_lat,
@@ -38,7 +38,6 @@ from summarize_nsrdb_county_polygon_archives import (
read_lookup, read_lookup,
) )
DEFAULT_ARCHIVE_DIR = Path("data/nrel/polygon_cloud_archives") DEFAULT_ARCHIVE_DIR = Path("data/nrel/polygon_cloud_archives")
DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json") DEFAULT_COUNTIES_GEOJSON = Path("data/geojson-counties-fips.json")
DEFAULT_REQUESTS_CSV = Path("data/nrel/county_polygon_cloud_request_manifest.csv") DEFAULT_REQUESTS_CSV = Path("data/nrel/county_polygon_cloud_request_manifest.csv")
@@ -59,15 +58,15 @@ SUMMARY_FIELDS = [
"rows_per_site_max", "rows_per_site_max",
"daylight_rows_per_site_min", "daylight_rows_per_site_min",
"daylight_rows_per_site_max", "daylight_rows_per_site_max",
"cloudinessIndexPct", "clearSkyGhiReductionIndex",
"areaWeightedCloudinessIndexPct", "areaWeightedClearSkyGhiReductionIndex",
"areaWeightedAvgObservedToClearskyRatio", "areaWeightedAvgObservedToClearskyRatio",
"areaWeightedSites", "areaWeightedSites",
"weightedCellAreaKm2", "weightedCellAreaKm2",
"countyAreaKm2", "countyAreaKm2",
"site_cloudiness_min", "site_clear_sky_ghi_reduction_min",
"site_cloudiness_max", "site_clear_sky_ghi_reduction_max",
"site_cloudiness_stddev", "site_clear_sky_ghi_reduction_stddev",
"clearOrProbablyClearPct", "clearOrProbablyClearPct",
"cloudyOrObscuredPct", "cloudyOrObscuredPct",
"fogPct", "fogPct",
@@ -75,7 +74,7 @@ SUMMARY_FIELDS = [
"iceCloudPct", "iceCloudPct",
"cirrusPct", "cirrusPct",
"unknownCloudTypePct", "unknownCloudTypePct",
"representativePointCloudinessIndexPct", "representativePointClearSkyGhiReductionIndex",
"areaWeightedMinusRepresentativePoint", "areaWeightedMinusRepresentativePoint",
"areaWeightedPctDiffFromRepresentativePoint", "areaWeightedPctDiffFromRepresentativePoint",
] ]
@@ -137,7 +136,7 @@ def pct(count: int, total: int) -> float:
def site_cloud_metrics(csv_text: str, min_clearsky_ghi: float) -> dict[str, float | int]: def site_cloud_metrics(csv_text: str, min_clearsky_ghi: float) -> dict[str, float | int]:
"""Return site-level cloudiness metrics from one NSRDB CSV response.""" """Return site-level clear-sky GHI reduction metrics from one NSRDB CSV response."""
rows = list(csv.reader(csv_text.splitlines())) rows = list(csv.reader(csv_text.splitlines()))
header_index = find_data_header(rows) header_index = find_data_header(rows)
if header_index is None: if header_index is None:
@@ -187,7 +186,7 @@ def site_cloud_metrics(csv_text: str, min_clearsky_ghi: float) -> dict[str, floa
raise ValueError("No daylight rows with valid GHI and Clearsky GHI were found.") raise ValueError("No daylight rows with valid GHI and Clearsky GHI were found.")
avg_ratio = ratio_sum / daylight_rows avg_ratio = ratio_sum / daylight_rows
cloudiness = 1 - avg_ratio clear_sky_ghi_reduction = 1 - avg_ratio
clear_or_probably_clear = daylight_cloud_type_counts.get("0", 0) + daylight_cloud_type_counts.get("1", 0) clear_or_probably_clear = daylight_cloud_type_counts.get("0", 0) + daylight_cloud_type_counts.get("1", 0)
cloudy_or_obscured = sum( cloudy_or_obscured = sum(
daylight_cloud_type_counts.get(code, 0) daylight_cloud_type_counts.get(code, 0)
@@ -206,7 +205,7 @@ def site_cloud_metrics(csv_text: str, min_clearsky_ghi: float) -> dict[str, floa
unknown_cloud_type = daylight_cloud_type_counts.get("10", 0) + daylight_cloud_type_counts.get("-15", 0) unknown_cloud_type = daylight_cloud_type_counts.get("10", 0) + daylight_cloud_type_counts.get("-15", 0)
return { return {
"cloudiness": cloudiness, "clear_sky_ghi_reduction": clear_sky_ghi_reduction,
"avg_ratio": avg_ratio, "avg_ratio": avg_ratio,
"all_rows": all_rows, "all_rows": all_rows,
"daylight_rows": daylight_rows, "daylight_rows": daylight_rows,
@@ -262,11 +261,11 @@ def summarize_archives(
if not site_metrics: if not site_metrics:
raise ValueError(f"No CSV files found in {', '.join(str(path) for path in paths)}") raise ValueError(f"No CSV files found in {', '.join(str(path) for path in paths)}")
cloudiness_values = [float(metrics["cloudiness"]) for metrics in site_metrics] reduction_values = [float(metrics["clear_sky_ghi_reduction"]) for metrics in site_metrics]
row_counts = [int(metrics["all_rows"]) for metrics in site_metrics] row_counts = [int(metrics["all_rows"]) for metrics in site_metrics]
daylight_row_counts = [int(metrics["daylight_rows"]) for metrics in site_metrics] daylight_row_counts = [int(metrics["daylight_rows"]) for metrics in site_metrics]
weighted_cloudiness = weighted_metric(site_metrics, "cloudiness", site_lon_lats, county_geometry, cell_size_m) weighted_reduction = weighted_metric(site_metrics, "clear_sky_ghi_reduction", site_lon_lats, county_geometry, cell_size_m)
weighted_avg_ratio = weighted_metric(site_metrics, "avg_ratio", site_lon_lats, county_geometry, cell_size_m) weighted_avg_ratio = weighted_metric(site_metrics, "avg_ratio", site_lon_lats, county_geometry, cell_size_m)
return { return {
@@ -275,14 +274,14 @@ def summarize_archives(
"rows_per_site_max": max(row_counts), "rows_per_site_max": max(row_counts),
"daylight_rows_per_site_min": min(daylight_row_counts), "daylight_rows_per_site_min": min(daylight_row_counts),
"daylight_rows_per_site_max": max(daylight_row_counts), "daylight_rows_per_site_max": max(daylight_row_counts),
"area_weighted_cloudiness": weighted_cloudiness["area_weighted_avg"], "area_weighted_clear_sky_ghi_reduction": weighted_reduction["area_weighted_avg"],
"area_weighted_avg_ratio": weighted_avg_ratio["area_weighted_avg"], "area_weighted_avg_ratio": weighted_avg_ratio["area_weighted_avg"],
"area_weighted_sites": weighted_cloudiness["area_weighted_sites"], "area_weighted_sites": weighted_reduction["area_weighted_sites"],
"weighted_cell_area_km2": weighted_cloudiness["weighted_cell_area_km2"], "weighted_cell_area_km2": weighted_reduction["weighted_cell_area_km2"],
"county_area_km2": weighted_cloudiness["county_area_km2"], "county_area_km2": weighted_reduction["county_area_km2"],
"site_cloudiness_min": min(cloudiness_values), "site_clear_sky_ghi_reduction_min": min(reduction_values),
"site_cloudiness_max": max(cloudiness_values), "site_clear_sky_ghi_reduction_max": max(reduction_values),
"site_cloudiness_stddev": statistics.pstdev(cloudiness_values) if len(cloudiness_values) > 1 else 0.0, "site_clear_sky_ghi_reduction_stddev": statistics.pstdev(reduction_values) if len(reduction_values) > 1 else 0.0,
"clear_or_probably_clear_pct": weighted_metric( "clear_or_probably_clear_pct": weighted_metric(
site_metrics, "clear_or_probably_clear_pct", site_lon_lats, county_geometry, cell_size_m site_metrics, "clear_or_probably_clear_pct", site_lon_lats, county_geometry, cell_size_m
)["area_weighted_avg"], )["area_weighted_avg"],
@@ -343,20 +342,20 @@ def build_summary_row(
if str(row.get("site_count", "")).strip().isdigit() if str(row.get("site_count", "")).strip().isdigit()
] ]
archive_zip = ";".join(str(path) for path in archive_paths) archive_zip = ";".join(str(path) for path in archive_paths)
area_weighted_cloudiness = float(polygon_summary["area_weighted_cloudiness"]) area_weighted_clear_sky_ghi_reduction = float(polygon_summary["area_weighted_clear_sky_ghi_reduction"])
representative_point_cloudiness = ( representative_point_reduction = (
float(representative_point_row["cloudinessIndexPct"]) float(representative_point_row["clearSkyGhiReductionIndex"])
if representative_point_row.get("cloudinessIndexPct") if representative_point_row.get("clearSkyGhiReductionIndex")
else None else None
) )
weighted_minus_representative_point = ( weighted_minus_representative_point = (
area_weighted_cloudiness - representative_point_cloudiness area_weighted_clear_sky_ghi_reduction - representative_point_reduction
if representative_point_cloudiness is not None if representative_point_reduction is not None
else None else None
) )
weighted_pct_diff_representative_point = ( weighted_pct_diff_representative_point = (
weighted_minus_representative_point / representative_point_cloudiness * 100 weighted_minus_representative_point / representative_point_reduction * 100
if representative_point_cloudiness not in (None, 0) if representative_point_reduction not in (None, 0)
else None else None
) )
@@ -371,15 +370,15 @@ def build_summary_row(
"rows_per_site_max": str(polygon_summary["rows_per_site_max"]), "rows_per_site_max": str(polygon_summary["rows_per_site_max"]),
"daylight_rows_per_site_min": str(polygon_summary["daylight_rows_per_site_min"]), "daylight_rows_per_site_min": str(polygon_summary["daylight_rows_per_site_min"]),
"daylight_rows_per_site_max": str(polygon_summary["daylight_rows_per_site_max"]), "daylight_rows_per_site_max": str(polygon_summary["daylight_rows_per_site_max"]),
"cloudinessIndexPct": format_float(area_weighted_cloudiness, 4), "clearSkyGhiReductionIndex": format_float(area_weighted_clear_sky_ghi_reduction, 4),
"areaWeightedCloudinessIndexPct": format_float(area_weighted_cloudiness, 4), "areaWeightedClearSkyGhiReductionIndex": format_float(area_weighted_clear_sky_ghi_reduction, 4),
"areaWeightedAvgObservedToClearskyRatio": format_float(float(polygon_summary["area_weighted_avg_ratio"]), 4), "areaWeightedAvgObservedToClearskyRatio": format_float(float(polygon_summary["area_weighted_avg_ratio"]), 4),
"areaWeightedSites": str(polygon_summary["area_weighted_sites"]), "areaWeightedSites": str(polygon_summary["area_weighted_sites"]),
"weightedCellAreaKm2": format_float(float(polygon_summary["weighted_cell_area_km2"]), 1), "weightedCellAreaKm2": format_float(float(polygon_summary["weighted_cell_area_km2"]), 1),
"countyAreaKm2": format_float(float(polygon_summary["county_area_km2"]), 1), "countyAreaKm2": format_float(float(polygon_summary["county_area_km2"]), 1),
"site_cloudiness_min": format_float(float(polygon_summary["site_cloudiness_min"]), 4), "site_clear_sky_ghi_reduction_min": format_float(float(polygon_summary["site_clear_sky_ghi_reduction_min"]), 4),
"site_cloudiness_max": format_float(float(polygon_summary["site_cloudiness_max"]), 4), "site_clear_sky_ghi_reduction_max": format_float(float(polygon_summary["site_clear_sky_ghi_reduction_max"]), 4),
"site_cloudiness_stddev": format_float(float(polygon_summary["site_cloudiness_stddev"]), 4), "site_clear_sky_ghi_reduction_stddev": format_float(float(polygon_summary["site_clear_sky_ghi_reduction_stddev"]), 4),
"clearOrProbablyClearPct": format_float(float(polygon_summary["clear_or_probably_clear_pct"]), 2), "clearOrProbablyClearPct": format_float(float(polygon_summary["clear_or_probably_clear_pct"]), 2),
"cloudyOrObscuredPct": format_float(float(polygon_summary["cloudy_or_obscured_pct"]), 2), "cloudyOrObscuredPct": format_float(float(polygon_summary["cloudy_or_obscured_pct"]), 2),
"fogPct": format_float(float(polygon_summary["fog_pct"]), 2), "fogPct": format_float(float(polygon_summary["fog_pct"]), 2),
@@ -387,7 +386,7 @@ def build_summary_row(
"iceCloudPct": format_float(float(polygon_summary["ice_cloud_pct"]), 2), "iceCloudPct": format_float(float(polygon_summary["ice_cloud_pct"]), 2),
"cirrusPct": format_float(float(polygon_summary["cirrus_pct"]), 2), "cirrusPct": format_float(float(polygon_summary["cirrus_pct"]), 2),
"unknownCloudTypePct": format_float(float(polygon_summary["unknown_cloud_type_pct"]), 2), "unknownCloudTypePct": format_float(float(polygon_summary["unknown_cloud_type_pct"]), 2),
"representativePointCloudinessIndexPct": format_float(representative_point_cloudiness, 4), "representativePointClearSkyGhiReductionIndex": format_float(representative_point_reduction, 4),
"areaWeightedMinusRepresentativePoint": format_float(weighted_minus_representative_point, 4), "areaWeightedMinusRepresentativePoint": format_float(weighted_minus_representative_point, 4),
"areaWeightedPctDiffFromRepresentativePoint": format_float(weighted_pct_diff_representative_point), "areaWeightedPctDiffFromRepresentativePoint": format_float(weighted_pct_diff_representative_point),
} }
@@ -465,8 +464,8 @@ def run(args: argparse.Namespace) -> None:
) )
print( print(
f"[{index}/{total}] {county_fips}: " f"[{index}/{total}] {county_fips}: "
f"sites={row['polygon_sites']}, weighted={row['areaWeightedCloudinessIndexPct']}, " f"sites={row['polygon_sites']}, weighted={row['areaWeightedClearSkyGhiReductionIndex']}, "
f"representative_point={row['representativePointCloudinessIndexPct'] or 'n/a'}" f"representative_point={row['representativePointClearSkyGhiReductionIndex'] or 'n/a'}"
f"{tile_note}{empty_note}" f"{tile_note}{empty_note}"
) )
@@ -500,7 +499,7 @@ def parse_args() -> argparse.Namespace:
"--min-clearsky-ghi", "--min-clearsky-ghi",
type=float, type=float,
default=DEFAULT_MIN_CLEARSKY_GHI, default=DEFAULT_MIN_CLEARSKY_GHI,
help="Minimum Clearsky GHI W/m2 for daylight cloudiness ratio rows.", help="Minimum Clearsky GHI W/m2 for daylight GHI ratio rows.",
) )
parser.add_argument( parser.add_argument(
"--reuse-existing-output", "--reuse-existing-output",
+97
View File
@@ -0,0 +1,97 @@
from __future__ import annotations
import csv
import sys
import unittest
from pathlib import Path
from tempfile import TemporaryDirectory
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR))
from apply_nsrdb_cloud_metric_to_climate_data import ( # noqa: E402
POLYGON_SOURCE_TAG,
REPRESENTATIVE_POINT_SOURCE_TAG,
merge_metric,
)
def write_csv(path: Path, fieldnames: list[str], rows: list[dict[str, str]]) -> None:
with path.open("w", encoding="utf-8", newline="") as handle:
writer = csv.DictWriter(handle, fieldnames=fieldnames)
writer.writeheader()
writer.writerows(rows)
class ApplyNsrdbCloudMetricTests(unittest.TestCase):
def test_polygon_area_weighted_value_wins_with_representative_fallback(self) -> None:
with TemporaryDirectory() as temp_dir:
base = Path(temp_dir)
climate_data = base / "climate-data.csv"
polygon_summary = base / "polygon-cloud.csv"
representative_summary = base / "representative-cloud.csv"
write_csv(
climate_data,
["countyFips", "meanDailyGlobalHorizontalRadiationKwhM2Day", "clearSkyGhiReductionIndex", "source"],
[
{
"countyFips": "01001",
"meanDailyGlobalHorizontalRadiationKwhM2Day": "4.8",
"clearSkyGhiReductionIndex": "0.9999",
"source": f"base + {REPRESENTATIVE_POINT_SOURCE_TAG}",
},
{
"countyFips": "01003",
"meanDailyGlobalHorizontalRadiationKwhM2Day": "4.9",
"clearSkyGhiReductionIndex": "",
"source": "base",
},
{
"countyFips": "01005",
"meanDailyGlobalHorizontalRadiationKwhM2Day": "5.0",
"clearSkyGhiReductionIndex": "",
"source": "base",
},
],
)
write_csv(
polygon_summary,
[
"county_fips",
"clearSkyGhiReductionIndex",
"areaWeightedClearSkyGhiReductionIndex",
],
[
{
"county_fips": "01001",
"clearSkyGhiReductionIndex": "0.1111",
"areaWeightedClearSkyGhiReductionIndex": "0.2222",
},
],
)
write_csv(
representative_summary,
["county_fips", "clearSkyGhiReductionIndex"],
[
{"county_fips": "01001", "clearSkyGhiReductionIndex": "0.3333"},
{"county_fips": "01003", "clearSkyGhiReductionIndex": "0.4444"},
],
)
result = merge_metric(climate_data, polygon_summary, representative_summary)
self.assertEqual(result, (3, 1, 1, 1))
with climate_data.open("r", encoding="utf-8", newline="") as handle:
rows = {row["countyFips"]: row for row in csv.DictReader(handle)}
self.assertEqual(rows["01001"]["clearSkyGhiReductionIndex"], "0.2222")
self.assertIn(POLYGON_SOURCE_TAG, rows["01001"]["source"])
self.assertNotIn(REPRESENTATIVE_POINT_SOURCE_TAG, rows["01001"]["source"])
self.assertEqual(rows["01003"]["clearSkyGhiReductionIndex"], "0.4444")
self.assertIn(REPRESENTATIVE_POINT_SOURCE_TAG, rows["01003"]["source"])
self.assertEqual(rows["01005"]["clearSkyGhiReductionIndex"], "")
if __name__ == "__main__":
unittest.main()
+211
View File
@@ -0,0 +1,211 @@
from __future__ import annotations
import csv
import json
import sys
import unittest
from pathlib import Path
from tempfile import TemporaryDirectory
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR))
from check_climate_data import ( # noqa: E402
CHECK_BLANKS,
CHECK_COLUMNS,
CHECK_COUNTIES,
CHECK_CROSS,
CHECK_FORMAT,
CHECK_GEOMETRY,
CHECK_SOURCES,
CHECK_VALUES,
DEFAULT_CLIMATE_CSV,
DEFAULT_COUNTIES_GEOJSON,
DEFAULT_METRIC_SOURCES,
EXPECTED_COLUMNS,
run_checks,
)
NOAA_GRID_COLUMNS = (
"avgTempF",
"avgDiurnalTempRangeF",
"annualPrecipIn",
"seasonalityIndex",
"wettestPrecipMonth",
"driestPrecipMonth",
"absoluteExtremeDays",
"avgSummerSpecificHumidityGKg",
"humidHeatDays",
"humidHeatSourceFips",
)
def make_row(fips: str, state: str, **overrides: str) -> dict[str, str]:
row = {
"countyFips": fips,
"countyName": "Test",
"state": state,
"koppenZone": "Cfa",
"koppenPrimaryClass": "",
"koppenSecondaryClass": "",
"avgTempF": "60.0",
"avgDiurnalTempRangeF": "20.00",
"annualPrecipIn": "40.0",
"seasonalityIndex": "20",
"wettestPrecipMonth": "May",
"driestPrecipMonth": "October",
"absoluteExtremeDays": "10.0",
"meanDailyGlobalHorizontalRadiationKwhM2Day": "4.5",
"clearSkyGhiReductionIndex": "0.25",
"avgSummerSpecificHumidityGKg": "12.0",
"humidHeatDays": "30.0",
"humidHeatSourceFips": fips,
"humidHeatFipsAdjustment": "",
"source": "test",
}
row.update(overrides)
return row
def valid_rows() -> list[dict[str, str]]:
alaska = make_row("02013", "AK", koppenZone="Dfc", **{column: "" for column in NOAA_GRID_COLUMNS})
return [make_row("01001", "AL"), alaska]
class CheckClimateDataTests(unittest.TestCase):
def setUp(self) -> None:
self._temp_dir = TemporaryDirectory()
self.base = Path(self._temp_dir.name)
self.csv_path = self.base / "climate-data.csv"
self.geojson_path = self.base / "counties.json"
self.sources_path = self.base / "metric_sources.json"
self.write_geojson(["01001", "02013"])
self.sources_path.write_text(json.dumps({"schemaVersion": 1, "metrics": {}}), encoding="utf-8")
def tearDown(self) -> None:
self._temp_dir.cleanup()
def write_csv(self, rows: list[dict[str, str]], fieldnames: list[str] | None = None, encoding: str = "utf-8") -> None:
with self.csv_path.open("w", encoding=encoding, newline="") as handle:
writer = csv.DictWriter(handle, fieldnames=fieldnames or list(EXPECTED_COLUMNS))
writer.writeheader()
writer.writerows(rows)
def write_geojson(self, fips_codes: list[str]) -> None:
features = [
{"type": "Feature", "properties": {"id": fips, "STATE": fips[:2]}, "geometry": None}
for fips in fips_codes
]
self.geojson_path.write_text(json.dumps({"type": "FeatureCollection", "features": features}), encoding="utf-8")
def run_report(self):
return run_checks(self.csv_path, self.geojson_path, self.sources_path)
def test_valid_file_passes_with_allowed_alaska_blanks(self) -> None:
self.write_csv(valid_rows())
report = self.run_report()
self.assertTrue(report.ok, report.checks)
def test_mixed_koppen_county_with_stripe_classes_passes(self) -> None:
rows = valid_rows()
rows[0].update(koppenZone="Mixed", koppenPrimaryClass="Csb", koppenSecondaryClass="Dsb")
self.write_csv(rows)
self.assertTrue(self.run_report().ok)
def test_mixed_county_without_stripe_classes_is_reported(self) -> None:
rows = valid_rows()
rows[0]["koppenZone"] = "Mixed"
self.write_csv(rows)
problems = self.run_report().checks[CHECK_CROSS]
self.assertEqual(problems, ["01001: Mixed Koppen county needs koppenPrimaryClass and koppenSecondaryClass"])
def test_stripe_classes_on_predominant_county_are_reported(self) -> None:
rows = valid_rows()
rows[0].update(koppenPrimaryClass="Csb", koppenSecondaryClass="Dsb")
self.write_csv(rows)
problems = self.run_report().checks[CHECK_CROSS]
self.assertEqual(problems, ["01001: Koppen stripe classes are set but koppenZone is Cfa, not Mixed"])
def test_invalid_or_identical_stripe_classes_are_reported(self) -> None:
rows = valid_rows()
rows[0].update(koppenZone="Mixed", koppenPrimaryClass="Csb", koppenSecondaryClass="Csb")
rows[1].update(koppenZone="Mixed", koppenPrimaryClass="Xyz", koppenSecondaryClass="Dfc")
self.write_csv(rows)
problems = self.run_report().checks[CHECK_CROSS]
self.assertEqual(
problems,
["01001: Koppen stripe classes are both Csb", "02013: invalid Koppen stripe classes 'Xyz'/'Dfc'"],
)
def test_bad_values_and_categories_are_reported(self) -> None:
rows = valid_rows()
rows[0].update(avgTempF="120", koppenZone="cfa", wettestPrecipMonth="Febuary", seasonalityIndex="12.5")
self.write_csv(rows)
problems = self.run_report().checks[CHECK_VALUES]
self.assertEqual(len(problems), 4, problems)
self.assertTrue(any("avgTempF=120 outside" in problem for problem in problems))
self.assertTrue(any("'cfa' is not an allowed category" in problem for problem in problems))
def test_blank_outside_allowed_states_is_reported(self) -> None:
rows = valid_rows()
rows[0]["avgTempF"] = ""
self.write_csv(rows)
problems = self.run_report().checks[CHECK_BLANKS]
self.assertEqual(problems, ["01001 (AL): avgTempF is blank"])
def test_county_missing_from_map_is_reported(self) -> None:
self.write_csv(valid_rows())
self.write_geojson(["01001"])
problems = self.run_report().checks[CHECK_GEOMETRY]
self.assertEqual(problems, ["02013 is in the CSV but has no map polygon"])
def test_duplicate_fips_and_unexpected_column_are_reported(self) -> None:
rows = valid_rows()
rows[1] = make_row("01001", "AL", extraMetric="1")
self.write_csv(rows, fieldnames=list(EXPECTED_COLUMNS) + ["extraMetric"])
report = self.run_report()
self.assertIn("duplicate countyFips 01001", report.checks[CHECK_COUNTIES])
self.assertTrue(any("unexpected column 'extraMetric'" in problem for problem in report.checks[CHECK_COLUMNS]))
def test_byte_order_mark_is_reported(self) -> None:
self.write_csv(valid_rows(), encoding="utf-8-sig")
report = self.run_report()
self.assertEqual(len(report.checks[CHECK_FORMAT]), 1)
self.assertEqual(report.checks[CHECK_COLUMNS], [])
def test_unknown_metric_in_sources_file_is_reported(self) -> None:
self.write_csv(valid_rows())
self.sources_path.write_text(json.dumps({"metrics": {"notAMetric": {}}}), encoding="utf-8")
problems = self.run_report().checks[CHECK_SOURCES]
self.assertEqual(problems, ["metric_sources.json describes unknown metric 'notAMetric'"])
@unittest.skipUnless(DEFAULT_CLIMATE_CSV.exists(), "project climate CSV not present")
def test_project_climate_csv_passes(self) -> None:
report = run_checks(DEFAULT_CLIMATE_CSV, DEFAULT_COUNTIES_GEOJSON, DEFAULT_METRIC_SOURCES)
self.assertTrue(report.ok, report.checks)
if __name__ == "__main__":
unittest.main()
+55
View File
@@ -0,0 +1,55 @@
from __future__ import annotations
import sys
import unittest
from pathlib import Path
import numpy as np
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR))
from summarize_county_gridmet_humidity import ( # noqa: E402
DEFAULT_HEAT_INDEX_THRESHOLD_F,
_heat_index_f,
)
class HeatIndexTests(unittest.TestCase):
def test_matches_nws_example(self) -> None:
result = _heat_index_f(np.array([100.0]), np.array([55.0]))
self.assertAlmostEqual(result[0], 124.0, delta=0.5)
def test_uses_simple_formula_below_regression_range(self) -> None:
result = _heat_index_f(np.array([70.0]), np.array([50.0]))
self.assertAlmostEqual(result[0], 69.525, places=3)
def test_default_threshold_is_extreme_caution_boundary(self) -> None:
heat_index = _heat_index_f(
np.array([85.0, 90.0]),
np.array([40.0, 55.0]),
)
self.assertEqual(DEFAULT_HEAT_INDEX_THRESHOLD_F, 90.0)
self.assertEqual(
(heat_index >= DEFAULT_HEAT_INDEX_THRESHOLD_F).tolist(),
[False, True],
)
def test_relative_humidity_is_clamped_to_physical_range(self) -> None:
result = _heat_index_f(
np.array([90.0, 90.0]),
np.array([-5.0, 105.0]),
)
expected = _heat_index_f(
np.array([90.0, 90.0]),
np.array([0.0, 100.0]),
)
np.testing.assert_allclose(result, expected)
if __name__ == "__main__":
unittest.main()
+104
View File
@@ -0,0 +1,104 @@
from __future__ import annotations
import sys
import unittest
from pathlib import Path
from tempfile import TemporaryDirectory
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import from_origin
from shapely.geometry import MultiPolygon, Polygon, box
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR))
from build_county_climate_data import _touched_raster_values, _zonal_majority_class # noqa: E402
from common.county_zonal_stats import geometry_window, split_at_antimeridian # noqa: E402
from common.koppen_legend import DEFAULT_KOPPEN_CODE_MAP # noqa: E402
def islands_on_both_sides() -> MultiPolygon:
"""A county shaped like Aleutians West: islands at +178..+180 and -180..-179."""
return MultiPolygon([box(178.2, 51.2, 179.8, 51.8), box(-179.8, 51.2, -179.2, 51.8)])
class SplitAtAntimeridianTests(unittest.TestCase):
def test_geometry_on_one_side_is_unchanged(self) -> None:
geometry = box(-100.0, 30.0, -90.0, 40.0)
pieces = split_at_antimeridian(geometry)
self.assertEqual(len(pieces), 1)
self.assertTrue(pieces[0].equals(geometry))
def test_islands_are_grouped_by_side(self) -> None:
geometry = MultiPolygon(
[box(172.0, 51.0, 173.0, 52.0), box(179.0, 51.0, 179.5, 52.0), box(-179.0, 51.0, -178.0, 52.0)]
)
pieces = split_at_antimeridian(geometry)
self.assertEqual(len(pieces), 2)
bounds = sorted(piece.bounds for piece in pieces)
self.assertEqual(bounds[0], (-179.0, 51.0, -178.0, 52.0))
self.assertEqual(bounds[1], (172.0, 51.0, 179.5, 52.0))
self.assertAlmostEqual(sum(piece.area for piece in pieces), geometry.area)
def test_single_outline_crossing_the_line_is_cut(self) -> None:
# Encoded in -180..180 coordinates, this 2-degree-wide ring looks 358 degrees wide.
crossing = Polygon([(179.0, 50.0), (-179.0, 50.0), (-179.0, 51.0), (179.0, 51.0)])
pieces = split_at_antimeridian(crossing)
bounds = sorted(piece.bounds for piece in pieces)
self.assertEqual(bounds, [(-180.0, 50.0, -179.0, 51.0), (179.0, 50.0, 180.0, 51.0)])
self.assertAlmostEqual(sum(piece.area for piece in pieces), 2.0)
class WindowedKoppenTests(unittest.TestCase):
def setUp(self) -> None:
self._temp_dir = TemporaryDirectory()
self.raster_path = Path(self._temp_dir.name) / "koppen.tif"
# Global 1-degree raster; row 38 covers 51..52 N.
data = np.zeros((180, 360), dtype=np.uint8)
data[38, 358] = 29 # ET at 178..179 E
data[38, 359] = 29 # ET at 179..180 E
data[38, 0] = 15 # Cfb at 180..179 W
with rasterio.open(
self.raster_path,
"w",
driver="GTiff",
height=180,
width=360,
count=1,
dtype="uint8",
crs="EPSG:4326",
transform=from_origin(-180.0, 90.0, 1.0, 1.0),
nodata=0,
) as destination:
destination.write(data, 1)
def tearDown(self) -> None:
self._temp_dir.cleanup()
def test_each_side_is_read_from_a_small_window(self) -> None:
with rasterio.open(self.raster_path) as source:
windows = [geometry_window(source, piece) for piece in split_at_antimeridian(islands_on_both_sides())]
values = _touched_raster_values(source, islands_on_both_sides(), split_antimeridian=True)
self.assertEqual(len(windows), 2)
self.assertTrue(all(window.width <= 5 for window in windows))
self.assertEqual(sorted(values[values != 0].tolist()), [15, 29, 29])
def test_majority_class_counts_islands_on_both_sides(self) -> None:
counties = gpd.GeoDataFrame(geometry=[islands_on_both_sides(), box(-100.0, 30.0, -99.0, 31.0)], crs="EPSG:4326")
classes = _zonal_majority_class(self.raster_path, counties, DEFAULT_KOPPEN_CODE_MAP)
self.assertEqual(classes, ["ET", "Cfa"])
if __name__ == "__main__":
unittest.main()
+281
View File
@@ -0,0 +1,281 @@
from __future__ import annotations
import csv
import math
import sys
import unittest
from pathlib import Path
from tempfile import TemporaryDirectory
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import from_origin
from shapely.geometry import MultiPolygon, box
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR))
from apply_koppen_metric_to_climate_data import apply_koppen_metric # noqa: E402
from build_county_koppen_metric import ( # noqa: E402
MIXED_CLASS,
build_koppen_records,
classify,
rank_class_shares,
)
from common.county_zonal_stats import area_weighted_class_weights # noqa: E402
from common.koppen_legend import DEFAULT_KOPPEN_CODE_MAP # noqa: E402
CFA, CFB, DFB, ET = 14, 15, 26, 29
def write_global_raster(path: Path, cells: dict[tuple[int, int], int]) -> None:
"""Write a 1-degree global raster; row r covers latitudes 89 - r to 90 - r."""
data = np.zeros((180, 360), dtype=np.uint8)
for (row, col), code in cells.items():
data[row, col] = code
with rasterio.open(
path,
"w",
driver="GTiff",
height=180,
width=360,
count=1,
dtype="uint8",
crs="EPSG:4326",
transform=from_origin(-180.0, 90.0, 1.0, 1.0),
nodata=0,
) as destination:
destination.write(data, 1)
def write_csv(path: Path, fieldnames: list[str], rows: list[dict[str, str]]) -> None:
with path.open("w", encoding="utf-8", newline="") as handle:
writer = csv.DictWriter(handle, fieldnames=fieldnames)
writer.writeheader()
writer.writerows(rows)
class ClassifyTests(unittest.TestCase):
def test_clear_majority_is_predominant(self) -> None:
self.assertEqual(classify([("Dfb", 0.60), ("Dfa", 0.30)]), "Dfb")
def test_exact_cutoffs_are_predominant(self) -> None:
self.assertEqual(classify([("Csa", 0.50), ("BSk", 0.45)]), "Csa")
def test_majority_with_close_runner_up_is_mixed(self) -> None:
self.assertEqual(classify([("Dfb", 0.501), ("Dfa", 0.499)]), MIXED_CLASS)
def test_no_majority_is_mixed(self) -> None:
self.assertEqual(classify([("BSh", 0.29), ("Csa", 0.24), ("Dsb", 0.23)]), MIXED_CLASS)
def test_single_class_is_predominant(self) -> None:
self.assertEqual(classify([("ET", 1.0)]), "ET")
def test_no_valid_cells_is_blank(self) -> None:
self.assertEqual(classify([]), "")
class RankClassSharesTests(unittest.TestCase):
def test_shares_are_ranked_and_ties_go_to_smaller_code(self) -> None:
ranked = rank_class_shares({DFB: 2.0, CFA: 2.0, ET: 1.0}, DEFAULT_KOPPEN_CODE_MAP)
self.assertEqual([name for name, _ in ranked], ["Cfa", "Dfb", "ET"])
self.assertAlmostEqual(ranked[0][1], 0.4)
self.assertAlmostEqual(sum(share for _, share in ranked), 1.0)
def test_unknown_code_raises(self) -> None:
with self.assertRaises(ValueError):
rank_class_shares({99: 1.0}, DEFAULT_KOPPEN_CODE_MAP)
def test_no_weights_give_no_classes(self) -> None:
self.assertEqual(rank_class_shares({}, DEFAULT_KOPPEN_CODE_MAP), [])
class KoppenRasterTestCase(unittest.TestCase):
def setUp(self) -> None:
self._temp_dir = TemporaryDirectory()
self.raster_path = Path(self._temp_dir.name) / "koppen.tif"
write_global_raster(
self.raster_path,
{
(59, 80): CFA, # 30..31 N, 100..99 W
(59, 81): CFB, # 30..31 N, 99..98 W
(49, 80): DFB, # 40..41 N, 100..99 W; 99..98 W is ocean
(29, 80): DFB, # 60..61 N
(79, 80): CFA, # 10..11 N
(38, 358): ET, # 51..52 N, 178..179 E
(38, 359): ET, # 51..52 N, 179..180 E
(38, 0): CFB, # 51..52 N, 180..179 W
},
)
def tearDown(self) -> None:
self._temp_dir.cleanup()
def weights(self, geometry) -> dict[int, float]:
with rasterio.open(self.raster_path) as source:
return area_weighted_class_weights(source, geometry)
class AreaWeightedClassWeightsTests(KoppenRasterTestCase):
def test_partial_cell_counts_by_fraction_inside(self) -> None:
weights = self.weights(box(-100.0, 30.0, -98.5, 31.0))
self.assertAlmostEqual(weights[CFB] / weights[CFA], 0.5)
def test_cells_are_weighted_by_latitude(self) -> None:
weights = self.weights(MultiPolygon([box(-100.0, 60.0, -99.0, 61.0), box(-100.0, 10.0, -99.0, 11.0)]))
expected = math.cos(math.radians(60.5)) / math.cos(math.radians(10.5))
self.assertAlmostEqual(weights[DFB] / weights[CFA], expected)
def test_ocean_cells_are_excluded(self) -> None:
weights = self.weights(box(-100.0, 40.0, -98.0, 41.0))
self.assertEqual(list(weights), [DFB])
def test_islands_on_both_sides_of_the_date_line_are_combined(self) -> None:
weights = self.weights(MultiPolygon([box(178.0, 51.0, 180.0, 52.0), box(-180.0, 51.0, -179.0, 52.0)]))
self.assertAlmostEqual(weights[ET] / weights[CFB], 2.0)
def test_no_valid_cells_return_no_weights(self) -> None:
self.assertEqual(self.weights(box(-50.0, 30.0, -49.0, 31.0)), {})
class BuildKoppenRecordsTests(KoppenRasterTestCase):
def test_counties_are_predominant_mixed_or_blank(self) -> None:
counties = gpd.GeoDataFrame(
{"county_fips": ["00001", "00002", "00003"], "county_name": ["A", "B", "C"], "state": ["AA"] * 3},
geometry=[box(-100.0, 40.0, -99.0, 41.0), box(-100.0, 30.0, -98.0, 31.0), box(-50.0, 30.0, -49.0, 31.0)],
crs="EPSG:4326",
)
records = build_koppen_records(counties, self.raster_path, DEFAULT_KOPPEN_CODE_MAP)
self.assertEqual([record["koppenZone"] for record in records], ["Dfb", MIXED_CLASS, ""])
self.assertEqual((records[0]["koppenTopShare"], records[0]["koppenSecondClass"]), ("1.0000", ""))
self.assertEqual(
(records[1]["koppenTopClass"], records[1]["koppenTopShare"], records[1]["koppenSecondShare"]),
("Cfa", "0.5000", "0.5000"),
)
self.assertEqual(records[2]["koppenTopShare"], "")
METRIC_TEST_FIELDS = ["countyFips", "koppenZone", "koppenTopClass", "koppenTopShare", "koppenSecondClass"]
class ApplyKoppenMetricTests(unittest.TestCase):
def setUp(self) -> None:
self._temp_dir = TemporaryDirectory()
base = Path(self._temp_dir.name)
self.climate_data = base / "climate-data.csv"
self.metric = base / "koppen.csv"
self.climate_fields = ["countyFips", "countyName", "state", "koppenZone", "avgTempF", "source"]
write_csv(
self.climate_data,
self.climate_fields,
[
{"countyFips": "01001", "countyName": "A", "state": "AL", "koppenZone": "Cfa", "avgTempF": "64.5", "source": "s"},
{"countyFips": "01003", "countyName": "B", "state": "AL", "koppenZone": "Dfb", "avgTempF": "50.1", "source": "s"},
],
)
self.write_metric(
[
{"countyFips": "01001", "koppenZone": "Mixed", "koppenTopClass": "Csb", "koppenTopShare": "0.4800", "koppenSecondClass": "Dsb"},
{"countyFips": "01003", "koppenZone": "Dfb", "koppenTopClass": "Dfb", "koppenTopShare": "0.9000", "koppenSecondClass": "Dfc"},
]
)
def tearDown(self) -> None:
self._temp_dir.cleanup()
def write_metric(self, rows: list[dict[str, str]]) -> None:
write_csv(self.metric, METRIC_TEST_FIELDS, rows)
def read_climate(self) -> tuple[list[str], list[dict[str, str]]]:
with self.climate_data.open("r", encoding="utf-8", newline="") as handle:
reader = csv.DictReader(handle)
return list(reader.fieldnames or []), list(reader)
def test_stripe_columns_are_added_after_koppen_zone(self) -> None:
changes, added_columns = apply_koppen_metric(self.climate_data, self.metric, self.climate_data)
fields, _ = self.read_climate()
self.assertEqual(added_columns, ["koppenPrimaryClass", "koppenSecondaryClass"])
self.assertEqual(
fields,
["countyFips", "countyName", "state", "koppenZone", "koppenPrimaryClass", "koppenSecondaryClass", "avgTempF", "source"],
)
self.assertEqual(
changes,
[
("01001", "koppenZone", "Cfa", "Mixed"),
("01001", "koppenPrimaryClass", "", "Csb"),
("01001", "koppenSecondaryClass", "", "Dsb"),
],
)
def test_only_mixed_counties_get_stripe_classes(self) -> None:
apply_koppen_metric(self.climate_data, self.metric, self.climate_data)
_, rows = self.read_climate()
self.assertEqual(
[(row["koppenZone"], row["koppenPrimaryClass"], row["koppenSecondaryClass"]) for row in rows],
[("Mixed", "Csb", "Dsb"), ("Dfb", "", "")],
)
self.assertEqual([row["avgTempF"] for row in rows], ["64.5", "50.1"])
def test_existing_stripe_columns_are_updated_in_place(self) -> None:
apply_koppen_metric(self.climate_data, self.metric, self.climate_data)
self.write_metric(
[
{"countyFips": "01001", "koppenZone": "Cfa", "koppenTopClass": "Cfa", "koppenTopShare": "0.7000", "koppenSecondClass": "Dfb"},
{"countyFips": "01003", "koppenZone": "Dfb", "koppenTopClass": "Dfb", "koppenTopShare": "0.9000", "koppenSecondClass": "Dfc"},
]
)
changes, added_columns = apply_koppen_metric(self.climate_data, self.metric, self.climate_data)
fields, rows = self.read_climate()
self.assertEqual(added_columns, [])
self.assertEqual(fields.count("koppenPrimaryClass"), 1)
self.assertEqual((rows[0]["koppenZone"], rows[0]["koppenPrimaryClass"], rows[0]["koppenSecondaryClass"]), ("Cfa", "", ""))
self.assertEqual(len(changes), 3)
def test_dry_run_writes_nothing(self) -> None:
before = self.climate_data.read_bytes()
changes, added_columns = apply_koppen_metric(self.climate_data, self.metric, self.climate_data, dry_run=True)
self.assertEqual(len(changes), 3)
self.assertEqual(len(added_columns), 2)
self.assertEqual(self.climate_data.read_bytes(), before)
def test_missing_county_raises_without_writing(self) -> None:
self.write_metric(
[{"countyFips": "01001", "koppenZone": "Mixed", "koppenTopClass": "Csb", "koppenTopShare": "0.4800", "koppenSecondClass": "Dsb"}]
)
before = self.climate_data.read_bytes()
with self.assertRaises(ValueError):
apply_koppen_metric(self.climate_data, self.metric, self.climate_data)
self.assertEqual(self.climate_data.read_bytes(), before)
def test_mixed_county_without_classes_raises_without_writing(self) -> None:
self.write_metric(
[
{"countyFips": "01001", "koppenZone": "Mixed", "koppenTopClass": "", "koppenTopShare": "", "koppenSecondClass": ""},
{"countyFips": "01003", "koppenZone": "Dfb", "koppenTopClass": "Dfb", "koppenTopShare": "0.9000", "koppenSecondClass": "Dfc"},
]
)
before = self.climate_data.read_bytes()
with self.assertRaises(ValueError):
apply_koppen_metric(self.climate_data, self.metric, self.climate_data)
self.assertEqual(self.climate_data.read_bytes(), before)
if __name__ == "__main__":
unittest.main()
@@ -7,7 +7,6 @@ from pathlib import Path
from tempfile import TemporaryDirectory from tempfile import TemporaryDirectory
from unittest.mock import patch from unittest.mock import patch
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts" SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR)) sys.path.insert(0, str(SCRIPTS_DIR))
@@ -6,7 +6,6 @@ from pathlib import Path
from tempfile import TemporaryDirectory from tempfile import TemporaryDirectory
from unittest.mock import patch from unittest.mock import patch
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts" SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR)) sys.path.insert(0, str(SCRIPTS_DIR))
@@ -1,28 +1,27 @@
from __future__ import annotations from __future__ import annotations
import csv
import sys import sys
import unittest import unittest
import csv
from datetime import datetime, timezone from datetime import datetime, timezone
from pathlib import Path from pathlib import Path
from tempfile import TemporaryDirectory from tempfile import TemporaryDirectory
from unittest.mock import call, patch from unittest.mock import call, patch
SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts" SCRIPTS_DIR = Path(__file__).resolve().parents[1] / "scripts"
sys.path.insert(0, str(SCRIPTS_DIR)) sys.path.insert(0, str(SCRIPTS_DIR))
from request_nsrdb_county_polygon_archives import ( # noqa: E402 from request_nsrdb_county_polygon_archives import ( # noqa: E402
EXCEPTIONAL_POLYGON_WAIT,
LARGE_POLYGON_WAIT,
MEDIUM_POLYGON_WAIT,
SMALL_POLYGON_WAIT,
ArchiveQueueMonitor, ArchiveQueueMonitor,
CountyRequestEvent, CountyRequestEvent,
CountyRequestState, CountyRequestState,
CountyRequestStateMachine, CountyRequestStateMachine,
EXCEPTIONAL_POLYGON_WAIT,
InvalidCountyRequestTransition, InvalidCountyRequestTransition,
LARGE_POLYGON_WAIT,
LocalQueueCapacityError, LocalQueueCapacityError,
MEDIUM_POLYGON_WAIT,
SMALL_POLYGON_WAIT,
) )