Initial Climate Mood Analysis project
This commit is contained in:
Justin Fisher
@@ -0,0 +1,690 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
Aggregate downloaded gridMET humidity/heat files to county-level metrics.
|
||||
|
||||
Expected input layout:
|
||||
data/gridmet/<year>/sph.nc
|
||||
data/gridmet/<year>/rmax.nc
|
||||
data/gridmet/<year>/rmin.nc
|
||||
|
||||
Default output:
|
||||
data/gridmet/county_gridmet_humidity.csv
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import csv
|
||||
import math
|
||||
from collections.abc import Sequence
|
||||
from dataclasses import dataclass
|
||||
from pathlib import Path
|
||||
|
||||
import geopandas as gpd
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import xarray as xr
|
||||
|
||||
|
||||
STATE_FIPS_TO_ABBR = {
|
||||
"01": "AL",
|
||||
"04": "AZ",
|
||||
"05": "AR",
|
||||
"06": "CA",
|
||||
"08": "CO",
|
||||
"09": "CT",
|
||||
"10": "DE",
|
||||
"11": "DC",
|
||||
"12": "FL",
|
||||
"13": "GA",
|
||||
"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",
|
||||
}
|
||||
CONUS_STATE_FIPS = set(STATE_FIPS_TO_ABBR)
|
||||
REQUIRED_VARIABLES = ("sph", "rmax", "rmin")
|
||||
SUMMER_MONTHS = {6, 7, 8}
|
||||
NOAA_REGION_CODE_TO_FIPS_OVERRIDES = {
|
||||
"18511": "11001",
|
||||
}
|
||||
NOAA_TMAX_SOURCE_FIPS_OVERRIDES = {
|
||||
"46113": (
|
||||
"46102",
|
||||
"Shannon County, SD changed name/code to Oglala Lakota County, SD effective 2015-05-01.",
|
||||
),
|
||||
"51515": (
|
||||
"51019",
|
||||
"Bedford independent city, VA changed to town status and was added to Bedford County effective 2013-07-01.",
|
||||
),
|
||||
"51678": (
|
||||
"51163",
|
||||
"Lexington independent city, VA uses surrounding Rockbridge County, VA as NOAA tmax proxy because NOAA county daily files do not include Lexington city.",
|
||||
),
|
||||
}
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class CountyGridMap:
|
||||
counties: gpd.GeoDataFrame
|
||||
pixel_indices: np.ndarray
|
||||
county_indices: np.ndarray
|
||||
weights: np.ndarray
|
||||
sorted_pixel_indices: np.ndarray
|
||||
sorted_county_indices: np.ndarray
|
||||
sorted_weights: np.ndarray
|
||||
group_starts: np.ndarray
|
||||
group_county_indices: np.ndarray
|
||||
lat_dim: str
|
||||
lon_dim: str
|
||||
|
||||
|
||||
@dataclass
|
||||
class YearData:
|
||||
datasets: list[xr.Dataset]
|
||||
arrays: dict[str, xr.DataArray]
|
||||
|
||||
|
||||
def _normalize_fips(value: object, width: int = 5) -> str:
|
||||
text = str(value).strip()
|
||||
digits = "".join(ch for ch in text if ch.isdigit())
|
||||
return digits.zfill(width)[-width:] if digits else ""
|
||||
|
||||
|
||||
def _load_state_crosswalk(path: Path) -> dict[str, str]:
|
||||
if not path.exists():
|
||||
print(
|
||||
f"State crosswalk not found at {path}. Falling back to the first two "
|
||||
"NOAA region-code digits as state FIPS."
|
||||
)
|
||||
return {}
|
||||
|
||||
rows = list(csv.reader(path.open("r", encoding="utf-8-sig", newline="")))
|
||||
if not rows:
|
||||
return {}
|
||||
|
||||
header = [cell.strip().lower() for cell in rows[0]]
|
||||
ncei_idx = next((idx for idx, cell in enumerate(header) if "ncei" in cell and "code" in cell), None)
|
||||
fips_idx = next((idx for idx, cell in enumerate(header) if "fips" in cell and "code" in cell), None)
|
||||
data_rows = rows[1:] if ncei_idx is not None and fips_idx is not None else rows
|
||||
|
||||
mapping: dict[str, str] = {}
|
||||
for row in data_rows:
|
||||
if ncei_idx is not None and fips_idx is not None:
|
||||
if len(row) <= max(ncei_idx, fips_idx):
|
||||
continue
|
||||
ncei = _normalize_fips(row[ncei_idx], 2)
|
||||
fips = _normalize_fips(row[fips_idx], 2)
|
||||
else:
|
||||
codes = [_normalize_fips(cell, 2) for cell in row]
|
||||
codes = [code for code in codes if code]
|
||||
if len(codes) < 2:
|
||||
continue
|
||||
ncei, fips = codes[0], codes[1]
|
||||
if ncei and fips:
|
||||
mapping[ncei] = fips
|
||||
return mapping
|
||||
|
||||
|
||||
def _county_fips_from_noaa_region(region_code: object, state_crosswalk: dict[str, str]) -> str:
|
||||
code = _normalize_fips(region_code, 5)
|
||||
if code in NOAA_REGION_CODE_TO_FIPS_OVERRIDES:
|
||||
return NOAA_REGION_CODE_TO_FIPS_OVERRIDES[code]
|
||||
state_fips = state_crosswalk.get(code[:2], code[:2])
|
||||
return f"{state_fips}{code[-3:]}" if state_fips else ""
|
||||
|
||||
|
||||
def _load_counties(path: Path) -> gpd.GeoDataFrame:
|
||||
counties = gpd.read_file(path)
|
||||
if counties.crs is None:
|
||||
counties = counties.set_crs("EPSG:4326")
|
||||
else:
|
||||
counties = counties.to_crs("EPSG:4326")
|
||||
|
||||
if "id" in counties.columns:
|
||||
fips = counties["id"]
|
||||
elif "GEOID" in counties.columns:
|
||||
fips = counties["GEOID"]
|
||||
elif "GEOID10" in counties.columns:
|
||||
fips = counties["GEOID10"]
|
||||
else:
|
||||
raise ValueError("County GeoJSON needs an id, GEOID, or GEOID10 column.")
|
||||
|
||||
counties["countyFips"] = fips.map(_normalize_fips)
|
||||
counties["stateFips"] = counties["countyFips"].str.slice(0, 2)
|
||||
counties = counties[counties["stateFips"].isin(CONUS_STATE_FIPS)].copy()
|
||||
counties["countyName"] = counties.get("NAME", counties["countyFips"]).astype(str)
|
||||
counties["state"] = counties["stateFips"].map(STATE_FIPS_TO_ABBR)
|
||||
counties = counties.sort_values("countyFips").reset_index(drop=True)
|
||||
counties["countyIndex"] = np.arange(len(counties), dtype=np.int32)
|
||||
return counties
|
||||
|
||||
|
||||
def _find_dimension(dataset: xr.Dataset, candidates: tuple[str, ...]) -> str:
|
||||
for candidate in candidates:
|
||||
if candidate in dataset.dims or candidate in dataset.coords:
|
||||
return candidate
|
||||
lower_lookup = {name.lower(): name for name in set(dataset.dims) | set(dataset.coords)}
|
||||
for candidate in candidates:
|
||||
if candidate.lower() in lower_lookup:
|
||||
return lower_lookup[candidate.lower()]
|
||||
raise ValueError(f"Unable to find one of these dimensions: {', '.join(candidates)}")
|
||||
|
||||
|
||||
def _select_variable(dataset: xr.Dataset, preferred: str) -> str:
|
||||
if preferred in dataset.data_vars:
|
||||
return preferred
|
||||
for name in dataset.data_vars:
|
||||
if name.lower() == preferred.lower():
|
||||
return name
|
||||
if len(dataset.data_vars) == 1:
|
||||
return next(iter(dataset.data_vars))
|
||||
raise ValueError(f"Unable to select variable {preferred}; found {list(dataset.data_vars)}")
|
||||
|
||||
|
||||
def _year_folder(path: Path) -> int | None:
|
||||
try:
|
||||
return int(path.name)
|
||||
except ValueError:
|
||||
return None
|
||||
|
||||
|
||||
def _discover_complete_years(gridmet_dir: Path) -> list[int]:
|
||||
years: list[int] = []
|
||||
for child in sorted(gridmet_dir.iterdir()):
|
||||
if not child.is_dir():
|
||||
continue
|
||||
year = _year_folder(child)
|
||||
if year is None:
|
||||
continue
|
||||
if all((child / f"{variable}.nc").exists() for variable in REQUIRED_VARIABLES):
|
||||
years.append(year)
|
||||
return years
|
||||
|
||||
|
||||
def _load_noaa_tmax_month(
|
||||
noaa_tmax_dir: Path,
|
||||
year: int,
|
||||
month: int,
|
||||
state_crosswalk: dict[str, str],
|
||||
) -> dict[str, list[float]]:
|
||||
path = noaa_tmax_dir / str(year) / f"tmax-{year}{month:02d}-cty-scaled.csv"
|
||||
if not path.exists():
|
||||
return {}
|
||||
|
||||
values: dict[str, list[float]] = {}
|
||||
with path.open("r", encoding="utf-8-sig", newline="") as handle:
|
||||
reader = csv.reader(handle)
|
||||
for row in reader:
|
||||
if len(row) < 7 or row[0].strip().lower() != "cty":
|
||||
continue
|
||||
county_fips = _county_fips_from_noaa_region(row[1], state_crosswalk)
|
||||
day_values: list[float] = []
|
||||
for raw_value in row[6:]:
|
||||
try:
|
||||
value = float(raw_value)
|
||||
except ValueError:
|
||||
day_values.append(math.nan)
|
||||
continue
|
||||
day_values.append(value if value > -999 else math.nan)
|
||||
values[county_fips] = day_values
|
||||
return values
|
||||
|
||||
|
||||
def _noaa_tmax_for_date(
|
||||
*,
|
||||
date: pd.Timestamp,
|
||||
county_fips: Sequence[str],
|
||||
noaa_tmax_dir: Path,
|
||||
state_crosswalk: dict[str, str],
|
||||
month_cache: dict[tuple[int, int], dict[str, list[float]]],
|
||||
) -> np.ndarray:
|
||||
cache_key = (int(date.year), int(date.month))
|
||||
if cache_key not in month_cache:
|
||||
month_cache[cache_key] = _load_noaa_tmax_month(
|
||||
noaa_tmax_dir,
|
||||
int(date.year),
|
||||
int(date.month),
|
||||
state_crosswalk,
|
||||
)
|
||||
|
||||
month_values = month_cache[cache_key]
|
||||
result = np.full(len(county_fips), np.nan, dtype=np.float64)
|
||||
day_offset = int(date.day) - 1
|
||||
for index, fips in enumerate(county_fips):
|
||||
source_fips = NOAA_TMAX_SOURCE_FIPS_OVERRIDES.get(fips, (fips, ""))[0]
|
||||
values = month_values.get(source_fips)
|
||||
if values is not None and day_offset < len(values):
|
||||
result[index] = values[day_offset]
|
||||
return result
|
||||
|
||||
|
||||
def _build_county_grid_map(counties: gpd.GeoDataFrame, sample_file: Path) -> CountyGridMap:
|
||||
with xr.open_dataset(sample_file) as dataset:
|
||||
lat_dim = _find_dimension(dataset, ("lat", "latitude", "y"))
|
||||
lon_dim = _find_dimension(dataset, ("lon", "longitude", "x"))
|
||||
lat_values = np.asarray(dataset[lat_dim].values)
|
||||
lon_values = np.asarray(dataset[lon_dim].values)
|
||||
|
||||
if lat_values.ndim != 1 or lon_values.ndim != 1:
|
||||
raise ValueError("This script expects 1D latitude and longitude coordinates.")
|
||||
|
||||
lon_grid, lat_grid = np.meshgrid(lon_values, lat_values)
|
||||
flat_lats = lat_grid.ravel()
|
||||
flat_lons = lon_grid.ravel()
|
||||
flat_indices = np.arange(flat_lats.size, dtype=np.int64)
|
||||
|
||||
bounds = counties.total_bounds
|
||||
in_bounds = (
|
||||
(flat_lons >= bounds[0] - 0.25)
|
||||
& (flat_lats >= bounds[1] - 0.25)
|
||||
& (flat_lons <= bounds[2] + 0.25)
|
||||
& (flat_lats <= bounds[3] + 0.25)
|
||||
)
|
||||
point_frame = gpd.GeoDataFrame(
|
||||
{
|
||||
"pixelIndex": flat_indices[in_bounds],
|
||||
"lat": flat_lats[in_bounds],
|
||||
},
|
||||
geometry=gpd.points_from_xy(flat_lons[in_bounds], flat_lats[in_bounds]),
|
||||
crs="EPSG:4326",
|
||||
)
|
||||
|
||||
joined = gpd.sjoin(
|
||||
point_frame,
|
||||
counties[["countyIndex", "geometry"]],
|
||||
how="inner",
|
||||
predicate="within",
|
||||
)
|
||||
joined = joined.drop_duplicates(subset=["pixelIndex", "countyIndex"])
|
||||
|
||||
pixel_indices = joined["pixelIndex"].to_numpy(dtype=np.int64)
|
||||
county_indices = joined["countyIndex"].to_numpy(dtype=np.int32)
|
||||
weights = np.cos(np.deg2rad(joined["lat"].to_numpy(dtype=np.float64)))
|
||||
|
||||
mapped_counties = set(county_indices.tolist())
|
||||
fallback_pixels: list[int] = []
|
||||
fallback_counties: list[int] = []
|
||||
fallback_weights: list[float] = []
|
||||
for county_index, geometry in enumerate(counties.geometry):
|
||||
if county_index in mapped_counties or geometry.is_empty:
|
||||
continue
|
||||
point = geometry.representative_point()
|
||||
lon_position = int(np.argmin(np.abs(lon_values - point.x)))
|
||||
lat_position = int(np.argmin(np.abs(lat_values - point.y)))
|
||||
flat_index = lat_position * len(lon_values) + lon_position
|
||||
fallback_pixels.append(flat_index)
|
||||
fallback_counties.append(county_index)
|
||||
fallback_weights.append(float(math.cos(math.radians(float(lat_values[lat_position])))))
|
||||
|
||||
if fallback_pixels:
|
||||
pixel_indices = np.concatenate([pixel_indices, np.asarray(fallback_pixels, dtype=np.int64)])
|
||||
county_indices = np.concatenate([county_indices, np.asarray(fallback_counties, dtype=np.int32)])
|
||||
weights = np.concatenate([weights, np.asarray(fallback_weights, dtype=np.float64)])
|
||||
|
||||
print(
|
||||
f"mapped {len(np.unique(county_indices)):,} counties to {len(pixel_indices):,} grid cells "
|
||||
f"({len(fallback_pixels):,} nearest-cell fallback counties)"
|
||||
)
|
||||
sort_order = np.argsort(county_indices, kind="stable")
|
||||
sorted_county_indices = county_indices[sort_order]
|
||||
group_starts = np.r_[0, np.flatnonzero(np.diff(sorted_county_indices)) + 1]
|
||||
group_county_indices = sorted_county_indices[group_starts]
|
||||
return CountyGridMap(
|
||||
counties=counties,
|
||||
pixel_indices=pixel_indices,
|
||||
county_indices=county_indices,
|
||||
weights=weights,
|
||||
sorted_pixel_indices=pixel_indices[sort_order],
|
||||
sorted_county_indices=sorted_county_indices,
|
||||
sorted_weights=weights[sort_order],
|
||||
group_starts=group_starts,
|
||||
group_county_indices=group_county_indices,
|
||||
lat_dim=lat_dim,
|
||||
lon_dim=lon_dim,
|
||||
)
|
||||
|
||||
|
||||
def _county_means(
|
||||
data_array: xr.DataArray,
|
||||
time_index: int,
|
||||
time_dim: str,
|
||||
county_map: CountyGridMap,
|
||||
) -> np.ndarray:
|
||||
data = data_array.isel({time_dim: time_index}).transpose(county_map.lat_dim, county_map.lon_dim)
|
||||
flat = np.asarray(data.values, dtype=np.float64).ravel()
|
||||
values = flat[county_map.pixel_indices]
|
||||
valid = np.isfinite(values)
|
||||
county_count = len(county_map.counties)
|
||||
sums = np.bincount(
|
||||
county_map.county_indices[valid],
|
||||
weights=values[valid] * county_map.weights[valid],
|
||||
minlength=county_count,
|
||||
)
|
||||
weight_sums = np.bincount(
|
||||
county_map.county_indices[valid],
|
||||
weights=county_map.weights[valid],
|
||||
minlength=county_count,
|
||||
)
|
||||
means = np.full(county_count, np.nan, dtype=np.float64)
|
||||
np.divide(sums, weight_sums, out=means, where=weight_sums > 0)
|
||||
return means
|
||||
|
||||
|
||||
def _county_means_chunk(
|
||||
data_array: xr.DataArray,
|
||||
start_index: int,
|
||||
end_index: int,
|
||||
time_dim: str,
|
||||
county_map: CountyGridMap,
|
||||
) -> np.ndarray:
|
||||
data = data_array.isel({time_dim: slice(start_index, end_index)}).transpose(
|
||||
time_dim,
|
||||
county_map.lat_dim,
|
||||
county_map.lon_dim,
|
||||
)
|
||||
raw = np.asarray(data.values, dtype=np.float64)
|
||||
selected = raw.reshape(raw.shape[0], -1)[:, county_map.sorted_pixel_indices]
|
||||
valid = np.isfinite(selected)
|
||||
|
||||
weighted_values = np.where(valid, selected * county_map.sorted_weights, 0.0)
|
||||
weight_values = np.where(valid, county_map.sorted_weights, 0.0)
|
||||
sums = np.add.reduceat(weighted_values, county_map.group_starts, axis=1)
|
||||
weight_sums = np.add.reduceat(weight_values, county_map.group_starts, axis=1)
|
||||
|
||||
grouped_means = np.full_like(sums, np.nan, dtype=np.float64)
|
||||
np.divide(sums, weight_sums, out=grouped_means, where=weight_sums > 0)
|
||||
|
||||
means = np.full((raw.shape[0], len(county_map.counties)), np.nan, dtype=np.float64)
|
||||
means[:, county_map.group_county_indices] = grouped_means
|
||||
return means
|
||||
|
||||
|
||||
def _heat_index_f(t_f: np.ndarray, rh_pct: np.ndarray) -> np.ndarray:
|
||||
rh = np.clip(rh_pct, 0, 100)
|
||||
heat_index = (
|
||||
-42.379
|
||||
+ 2.04901523 * t_f
|
||||
+ 10.14333127 * rh
|
||||
- 0.22475541 * t_f * rh
|
||||
- 0.00683783 * t_f * t_f
|
||||
- 0.05481717 * rh * rh
|
||||
+ 0.00122874 * t_f * t_f * rh
|
||||
+ 0.00085282 * t_f * rh * rh
|
||||
- 0.00000199 * t_f * t_f * rh * rh
|
||||
)
|
||||
|
||||
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)
|
||||
heat_index = np.where((rh < 13) & (80 <= t_f) & (t_f <= 112), heat_index - low_rh_adjustment, heat_index)
|
||||
heat_index = np.where((rh > 85) & (80 <= t_f) & (t_f <= 87), heat_index + high_rh_adjustment, heat_index)
|
||||
return np.where(t_f >= 80, heat_index, t_f)
|
||||
|
||||
|
||||
def _time_dimension(data_array: xr.DataArray, lat_dim: str, lon_dim: str) -> str:
|
||||
candidates = [dim for dim in data_array.dims if dim not in {lat_dim, lon_dim}]
|
||||
if not candidates:
|
||||
raise ValueError(f"No time dimension found for {data_array.name}.")
|
||||
return candidates[0]
|
||||
|
||||
|
||||
def _open_year_datasets(gridmet_dir: Path, year: int) -> YearData:
|
||||
datasets: list[xr.Dataset] = []
|
||||
arrays: dict[str, xr.DataArray] = {}
|
||||
for variable in REQUIRED_VARIABLES:
|
||||
dataset = xr.open_dataset(gridmet_dir / str(year) / f"{variable}.nc")
|
||||
datasets.append(dataset)
|
||||
variable_name = _select_variable(dataset, variable)
|
||||
arrays[variable] = dataset[variable_name]
|
||||
return YearData(datasets=datasets, arrays=arrays)
|
||||
|
||||
|
||||
def _close_year_datasets(year_data: YearData) -> None:
|
||||
for dataset in year_data.datasets:
|
||||
try:
|
||||
dataset.close()
|
||||
except Exception:
|
||||
pass
|
||||
|
||||
|
||||
def _round_or_blank(value: float, digits: int) -> str:
|
||||
if not np.isfinite(value):
|
||||
return ""
|
||||
return str(round(float(value), digits))
|
||||
|
||||
|
||||
def summarize(args: argparse.Namespace) -> None:
|
||||
years = args.years or _discover_complete_years(args.gridmet_dir)
|
||||
if not years:
|
||||
raise FileNotFoundError(
|
||||
f"No complete gridMET years found in {args.gridmet_dir}. "
|
||||
"Run scripts/download_gridmet_data.py first."
|
||||
)
|
||||
|
||||
if args.years is None:
|
||||
start = args.start_year if args.start_year is not None else 1991
|
||||
end = args.end_year if args.end_year is not None else 2020
|
||||
years = [year for year in years if start <= year <= end]
|
||||
elif args.start_year is not None or args.end_year is not None:
|
||||
start = args.start_year if args.start_year is not None else min(years)
|
||||
end = args.end_year if args.end_year is not None else max(years)
|
||||
years = [year for year in years if start <= year <= end]
|
||||
|
||||
if not years:
|
||||
raise FileNotFoundError("No complete gridMET years matched the requested year filters.")
|
||||
|
||||
sample_file = args.gridmet_dir / str(years[0]) / "sph.nc"
|
||||
counties = _load_counties(args.counties_geojson)
|
||||
county_map = _build_county_grid_map(counties, sample_file)
|
||||
county_count = len(county_map.counties)
|
||||
county_fips = county_map.counties["countyFips"].astype(str).tolist()
|
||||
state_crosswalk = _load_state_crosswalk(args.state_crosswalk)
|
||||
noaa_month_cache: dict[tuple[int, int], dict[str, list[float]]] = {}
|
||||
|
||||
annual_sph_sum = np.zeros(county_count, dtype=np.float64)
|
||||
annual_sph_count = np.zeros(county_count, dtype=np.float64)
|
||||
summer_sph_sum = np.zeros(county_count, dtype=np.float64)
|
||||
summer_sph_count = np.zeros(county_count, dtype=np.float64)
|
||||
annual_rh_sum = np.zeros(county_count, dtype=np.float64)
|
||||
annual_rh_count = np.zeros(county_count, dtype=np.float64)
|
||||
summer_rh_sum = np.zeros(county_count, dtype=np.float64)
|
||||
summer_rh_count = np.zeros(county_count, dtype=np.float64)
|
||||
heat_day_sum = np.zeros(county_count, dtype=np.float64)
|
||||
years_with_heat_data = np.zeros(county_count, dtype=np.float64)
|
||||
|
||||
for year in years:
|
||||
print(f"process {year}")
|
||||
year_data = _open_year_datasets(args.gridmet_dir, year)
|
||||
arrays = year_data.arrays
|
||||
try:
|
||||
time_dim = _time_dimension(arrays["sph"], county_map.lat_dim, county_map.lon_dim)
|
||||
dates = pd.to_datetime(arrays["sph"][time_dim].values)
|
||||
yearly_heat_days = np.zeros(county_count, dtype=np.float64)
|
||||
yearly_heat_valid_days = np.zeros(county_count, dtype=np.float64)
|
||||
|
||||
for chunk_start in range(0, len(dates), args.chunk_days):
|
||||
chunk_end = min(chunk_start + args.chunk_days, len(dates))
|
||||
chunk_dates = dates[chunk_start:chunk_end]
|
||||
sph_chunk = _county_means_chunk(arrays["sph"], chunk_start, chunk_end, time_dim, county_map) * 1000
|
||||
rmax_chunk = _county_means_chunk(arrays["rmax"], chunk_start, chunk_end, time_dim, county_map)
|
||||
rmin_chunk = _county_means_chunk(arrays["rmin"], chunk_start, chunk_end, time_dim, county_map)
|
||||
|
||||
for offset, date in enumerate(chunk_dates):
|
||||
sph = sph_chunk[offset]
|
||||
rmax = rmax_chunk[offset]
|
||||
rmin = rmin_chunk[offset]
|
||||
tmax_f = _noaa_tmax_for_date(
|
||||
date=date,
|
||||
county_fips=county_fips,
|
||||
noaa_tmax_dir=args.noaa_tmax_dir,
|
||||
state_crosswalk=state_crosswalk,
|
||||
month_cache=noaa_month_cache,
|
||||
) * 9 / 5 + 32
|
||||
|
||||
rh_mean = (rmax + rmin) / 2
|
||||
annual_sph_valid = np.isfinite(sph)
|
||||
annual_rh_valid = np.isfinite(rh_mean)
|
||||
annual_sph_sum[annual_sph_valid] += sph[annual_sph_valid]
|
||||
annual_sph_count[annual_sph_valid] += 1
|
||||
annual_rh_sum[annual_rh_valid] += rh_mean[annual_rh_valid]
|
||||
annual_rh_count[annual_rh_valid] += 1
|
||||
|
||||
if int(date.month) in SUMMER_MONTHS:
|
||||
summer_sph_sum[annual_sph_valid] += sph[annual_sph_valid]
|
||||
summer_sph_count[annual_sph_valid] += 1
|
||||
summer_rh_sum[annual_rh_valid] += rh_mean[annual_rh_valid]
|
||||
summer_rh_count[annual_rh_valid] += 1
|
||||
|
||||
heat_valid = np.isfinite(tmax_f) & np.isfinite(rmin)
|
||||
heat_index = _heat_index_f(tmax_f, rmin)
|
||||
yearly_heat_days[heat_valid] += heat_index[heat_valid] >= args.heat_index_threshold_f
|
||||
yearly_heat_valid_days[heat_valid] += 1
|
||||
|
||||
has_heat_year = yearly_heat_valid_days > 0
|
||||
heat_day_sum[has_heat_year] += yearly_heat_days[has_heat_year]
|
||||
years_with_heat_data[has_heat_year] += 1
|
||||
finally:
|
||||
_close_year_datasets(year_data)
|
||||
|
||||
avg_sph = np.divide(annual_sph_sum, annual_sph_count, out=np.full(county_count, np.nan), where=annual_sph_count > 0)
|
||||
avg_summer_sph = np.divide(
|
||||
summer_sph_sum,
|
||||
summer_sph_count,
|
||||
out=np.full(county_count, np.nan),
|
||||
where=summer_sph_count > 0,
|
||||
)
|
||||
avg_rh = np.divide(annual_rh_sum, annual_rh_count, out=np.full(county_count, np.nan), where=annual_rh_count > 0)
|
||||
avg_summer_rh = np.divide(
|
||||
summer_rh_sum,
|
||||
summer_rh_count,
|
||||
out=np.full(county_count, np.nan),
|
||||
where=summer_rh_count > 0,
|
||||
)
|
||||
humid_heat_days = np.divide(
|
||||
heat_day_sum,
|
||||
years_with_heat_data,
|
||||
out=np.full(county_count, np.nan),
|
||||
where=years_with_heat_data > 0,
|
||||
)
|
||||
source_period = f"{min(years)}-{max(years)}" if len(years) > 1 else str(years[0])
|
||||
|
||||
args.output.parent.mkdir(parents=True, exist_ok=True)
|
||||
with args.output.open("w", newline="", encoding="utf-8") as handle:
|
||||
writer = csv.DictWriter(
|
||||
handle,
|
||||
fieldnames=[
|
||||
"countyFips",
|
||||
"countyName",
|
||||
"state",
|
||||
"avgSpecificHumidityGKg",
|
||||
"avgSummerSpecificHumidityGKg",
|
||||
"avgRelHumidityPct",
|
||||
"avgSummerRelHumidityPct",
|
||||
"humidHeatDays",
|
||||
"humidHeatSourceFips",
|
||||
"humidHeatFipsAdjustment",
|
||||
"gridmetAnalysisYears",
|
||||
"gridmetSource",
|
||||
],
|
||||
)
|
||||
writer.writeheader()
|
||||
for index, county in county_map.counties.iterrows():
|
||||
writer.writerow(
|
||||
{
|
||||
"countyFips": county["countyFips"],
|
||||
"countyName": county["countyName"],
|
||||
"state": county["state"],
|
||||
"avgSpecificHumidityGKg": _round_or_blank(avg_sph[index], 3),
|
||||
"avgSummerSpecificHumidityGKg": _round_or_blank(avg_summer_sph[index], 3),
|
||||
"avgRelHumidityPct": _round_or_blank(avg_rh[index], 2),
|
||||
"avgSummerRelHumidityPct": _round_or_blank(avg_summer_rh[index], 2),
|
||||
"humidHeatDays": _round_or_blank(humid_heat_days[index], 2),
|
||||
"humidHeatSourceFips": NOAA_TMAX_SOURCE_FIPS_OVERRIDES.get(
|
||||
county["countyFips"],
|
||||
(county["countyFips"], ""),
|
||||
)[0],
|
||||
"humidHeatFipsAdjustment": NOAA_TMAX_SOURCE_FIPS_OVERRIDES.get(
|
||||
county["countyFips"],
|
||||
("", ""),
|
||||
)[1],
|
||||
"gridmetAnalysisYears": len(years),
|
||||
"gridmetSource": (
|
||||
f"gridmet-{source_period} county-cell-weighted; "
|
||||
"relative-humidity=(rmax+rmin)/2; "
|
||||
"humidHeatDays=noaa-nclimgrid-tmax+gridmet-rmin heat-index-gte-90f"
|
||||
),
|
||||
}
|
||||
)
|
||||
|
||||
print(f"wrote {args.output}")
|
||||
|
||||
|
||||
def _parse_years(text: str) -> list[int]:
|
||||
years: set[int] = set()
|
||||
for part in text.split(","):
|
||||
token = part.strip()
|
||||
if not token:
|
||||
continue
|
||||
if "-" in token:
|
||||
start_text, end_text = token.split("-", 1)
|
||||
years.update(range(int(start_text), int(end_text) + 1))
|
||||
else:
|
||||
years.add(int(token))
|
||||
return sorted(years)
|
||||
|
||||
|
||||
def main() -> int:
|
||||
parser = argparse.ArgumentParser(
|
||||
description="Summarize downloaded gridMET humidity files into county-level CSV metrics."
|
||||
)
|
||||
parser.add_argument("--gridmet-dir", type=Path, default=Path("data/gridmet"))
|
||||
parser.add_argument("--noaa-tmax-dir", type=Path, default=Path("data/noaa/tmax_cty_scaled"))
|
||||
parser.add_argument("--state-crosswalk", type=Path, default=Path("data/noaa/nclimgrid/us-state-codes_ncei-to-fips.csv"))
|
||||
parser.add_argument("--counties-geojson", type=Path, default=Path("data/geojson-counties-fips.json"))
|
||||
parser.add_argument("--output", type=Path, default=Path("data/gridmet/county_gridmet_humidity.csv"))
|
||||
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("--end-year", type=int, default=None)
|
||||
parser.add_argument("--heat-index-threshold-f", type=float, default=90.0)
|
||||
parser.add_argument("--chunk-days", type=int, default=31)
|
||||
summarize(parser.parse_args())
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main())
|
||||
Reference in New Issue
Block a user