{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91ba68f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Download and unzip all 5-minute (05m) precipitation frequency files\n",
    "for Texas from NOAA HDSC Atlas 14.\n",
    "\n",
    "Files downloaded:\n",
    "  tx{N}yr05ma.zip   - point estimates\n",
    "  tx{N}yr05mal.zip  - lower confidence bound\n",
    "  tx{N}yr05mau.zip  - upper confidence bound\n",
    "\n",
    "for return periods: 1, 2, 5, 10, 25, 50, 100, 200, 500, 1000 years\n",
    "\"\"\"\n",
    "\n",
    "import os\n",
    "import urllib.request\n",
    "import zipfile\n",
    "\n",
    "# --- Configuration ---\n",
    "BASE_URL = \"https://hdsc.nws.noaa.gov/pub/hdsc/data/tx/\"\n",
    "OUTPUT_DIR = r\"E:\\texas_precip\"\n",
    "RETURN_PERIODS = [1, 2, 5, 10, 25, 50, 100, 200, 500, 1000]\n",
    "SUFFIXES = [\"a\"]  # point estimate only\n",
    "\n",
    "os.makedirs(OUTPUT_DIR, exist_ok=True)\n",
    "\n",
    "\n",
    "def download_and_unzip(url, dest_dir):\n",
    "    filename = url.split(\"/\")[-1]\n",
    "    filepath = os.path.join(dest_dir, filename)\n",
    "\n",
    "    print(f\"Downloading {filename} ...\", end=\" \", flush=True)\n",
    "    try:\n",
    "        urllib.request.urlretrieve(url, filepath)\n",
    "        print(\"OK\", end=\" \", flush=True)\n",
    "    except Exception as e:\n",
    "        print(f\"FAILED ({e})\")\n",
    "        return\n",
    "\n",
    "    try:\n",
    "        with zipfile.ZipFile(filepath, \"r\") as z:\n",
    "            z.extractall(dest_dir)\n",
    "        os.remove(filepath)  # remove zip after extracting\n",
    "        print(\"-> unzipped\")\n",
    "    except zipfile.BadZipFile:\n",
    "        print(\"-> not a valid zip, kept as-is\")\n",
    "\n",
    "\n",
    "for rp in RETURN_PERIODS:\n",
    "    for suffix in SUFFIXES:\n",
    "        # 1-year files use a slightly different naming pattern (no _ams variant here)\n",
    "        filename = f\"tx{rp}yr05m{suffix}.zip\"\n",
    "        url = BASE_URL + filename\n",
    "        download_and_unzip(url, OUTPUT_DIR)\n",
    "\n",
    "print(f\"\\nDone. Files saved to: {os.path.abspath(OUTPUT_DIR)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60124678",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Convert NOAA Atlas 14 Texas 5-minute ASC precipitation files\n",
    "to a single NetCDF4 file.\n",
    "\n",
    "- Input : all tx{N}yr05ma.asc files in INPUT_DIR\n",
    "- Output: tx_5min_intensity.nc in OUTPUT_DIR\n",
    "- Structure: single variable intensity(return_period, lat, lon)\n",
    "- Values: converted from integer (thousandths of inches depth)\n",
    "          to intensity in inches/hour (depth_in * 12)\n",
    "- Missing: -9 (raw) -> NaN\n",
    "\"\"\"\n",
    "\n",
    "import os\n",
    "import re\n",
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "\n",
    "# --- Configuration ---\n",
    "INPUT_DIR = r\"E:\\texas_precip\"\n",
    "OUTPUT_NC = r\"E:\\texas_precip\\tx_5min_intensity.nc\"\n",
    "\n",
    "RETURN_PERIODS = [1, 2, 5, 10, 25, 50, 100, 200, 500, 1000]\n",
    "NODATA_RAW   = -9\n",
    "SCALE        = 1000.0   # values stored as thousandths of inches\n",
    "DURATION_MIN = 5\n",
    "IN_HR_FACTOR = 60.0 / DURATION_MIN  # = 12.0\n",
    "\n",
    "\n",
    "def parse_asc_header(f):\n",
    "    header = {}\n",
    "    for _ in range(6):\n",
    "        line = f.readline().strip().split()\n",
    "        header[line[0].lower()] = float(line[1])\n",
    "    return header\n",
    "\n",
    "\n",
    "def read_asc(filepath):\n",
    "    with open(filepath, \"r\") as f:\n",
    "        hdr = parse_asc_header(f)\n",
    "        data = np.loadtxt(f, dtype=np.float32)\n",
    "    return hdr, data\n",
    "\n",
    "\n",
    "# --- Discover files ---\n",
    "pattern = re.compile(r\"tx(\\d+)yr05ma\\.asc$\", re.IGNORECASE)\n",
    "found = {}\n",
    "for fname in os.listdir(INPUT_DIR):\n",
    "    m = pattern.match(fname)\n",
    "    if m:\n",
    "        rp = int(m.group(1))\n",
    "        found[rp] = os.path.join(INPUT_DIR, fname)\n",
    "\n",
    "missing = [rp for rp in RETURN_PERIODS if rp not in found]\n",
    "if missing:\n",
    "    print(f\"WARNING: No ASC file found for return periods: {missing}\")\n",
    "\n",
    "rps_to_use = sorted([rp for rp in RETURN_PERIODS if rp in found])\n",
    "if not rps_to_use:\n",
    "    raise FileNotFoundError(f\"No matching ASC files found in {INPUT_DIR}\")\n",
    "\n",
    "print(f\"Found {len(rps_to_use)} ASC files: {rps_to_use}\")\n",
    "\n",
    "# --- Read first file to get grid geometry ---\n",
    "hdr0, _ = read_asc(found[rps_to_use[0]])\n",
    "ncols  = int(hdr0[\"ncols\"])\n",
    "nrows  = int(hdr0[\"nrows\"])\n",
    "xll    = hdr0[\"xllcorner\"]\n",
    "yll    = hdr0[\"yllcorner\"]\n",
    "cell   = hdr0[\"cellsize\"]\n",
    "nodata = hdr0.get(\"nodata_value\", NODATA_RAW)\n",
    "\n",
    "# Cell-center coordinates, lat north-up\n",
    "lons = xll + (np.arange(ncols) + 0.5) * cell\n",
    "lats = yll + (np.arange(nrows - 1, -1, -1) + 0.5) * cell\n",
    "\n",
    "# --- Write NetCDF ---\n",
    "print(f\"Writing {OUTPUT_NC} ...\")\n",
    "ds = nc.Dataset(OUTPUT_NC, \"w\", format=\"NETCDF4\")\n",
    "\n",
    "ds.title       = \"NOAA Atlas 14 Texas 5-Minute Precipitation Intensity\"\n",
    "ds.institution = \"NOAA/NWS Hydrometeorological Design Studies Center\"\n",
    "ds.source      = \"NOAA Atlas 14 Volume 11\"\n",
    "ds.Conventions = \"CF-1.8\"\n",
    "ds.crs         = \"NAD83 geographic (EPSG:4269)\"\n",
    "\n",
    "# Dimensions\n",
    "ds.createDimension(\"return_period\", len(rps_to_use))\n",
    "ds.createDimension(\"lat\", nrows)\n",
    "ds.createDimension(\"lon\", ncols)\n",
    "\n",
    "# Coordinate variables\n",
    "rp_var = ds.createVariable(\"return_period\", \"i4\", (\"return_period\",))\n",
    "rp_var.units     = \"years\"\n",
    "rp_var.long_name = \"return period\"\n",
    "rp_var[:]        = np.array(rps_to_use, dtype=np.int32)\n",
    "\n",
    "lat_var = ds.createVariable(\"lat\", \"f4\", (\"lat\",))\n",
    "lat_var.units         = \"degrees_north\"\n",
    "lat_var.standard_name = \"latitude\"\n",
    "lat_var[:]            = lats\n",
    "\n",
    "lon_var = ds.createVariable(\"lon\", \"f4\", (\"lon\",))\n",
    "lon_var.units         = \"degrees_east\"\n",
    "lon_var.standard_name = \"longitude\"\n",
    "lon_var[:]            = lons\n",
    "\n",
    "# Single intensity variable\n",
    "intensity_var = ds.createVariable(\n",
    "    \"intensity\", \"f4\", (\"return_period\", \"lat\", \"lon\"),\n",
    "    fill_value=np.nan, zlib=True, complevel=4\n",
    ")\n",
    "intensity_var.long_name     = \"5-minute precipitation intensity\"\n",
    "intensity_var.units         = \"inches/hour\"\n",
    "intensity_var.duration_min  = DURATION_MIN\n",
    "intensity_var.grid_mapping  = \"crs\"\n",
    "intensity_var.coordinates   = \"return_period lat lon\"\n",
    "\n",
    "# Fill intensity array\n",
    "for i, rp in enumerate(rps_to_use):\n",
    "    print(f\"  Processing {rp}-year ...\", end=\" \", flush=True)\n",
    "    hdr, raw = read_asc(found[rp])\n",
    "    data = raw.astype(np.float32)\n",
    "    data[data == nodata]     = np.nan\n",
    "    data[data == NODATA_RAW] = np.nan\n",
    "    intensity_var[i, :, :]  = (data / SCALE) * IN_HR_FACTOR\n",
    "    print(\"OK\")\n",
    "\n",
    "# CRS variable\n",
    "crs_var = ds.createVariable(\"crs\", \"i4\")\n",
    "crs_var.grid_mapping_name           = \"latitude_longitude\"\n",
    "crs_var.longitude_of_prime_meridian = 0.0\n",
    "crs_var.semi_major_axis             = 6378137.0\n",
    "crs_var.inverse_flattening          = 298.257222\n",
    "crs_var.crs_wkt = (\n",
    "    'GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",'\n",
    "    'SPHEROID[\"GRS 1980\",6378137,298.257222101]],'\n",
    "    'PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]'\n",
    ")\n",
    "\n",
    "ds.close()\n",
    "print(f\"\\nDone. Output: {OUTPUT_NC}\")\n",
    "print(f\"  Grid  : {nrows} rows x {ncols} cols\")\n",
    "print(f\"  Layers: {rps_to_use} year return periods\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef8ed4d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Download and unzip all 1-year precipitation frequency files\n",
    "for Texas from NOAA HDSC Atlas 14, up to 24-hour duration.\n",
    "\n",
    "Durations included:\n",
    "  5, 10, 15, 30, 45, 60 minutes\n",
    "  2, 3, 6, 12, 24 hours\n",
    "\n",
    "Only point estimates (*a.zip) are downloaded.\n",
    "Files saved to E:\\texas_precip\n",
    "\"\"\"\n",
    "\n",
    "import os\n",
    "import urllib.request\n",
    "import zipfile\n",
    "\n",
    "BASE_URL   = \"https://hdsc.nws.noaa.gov/pub/hdsc/data/tx/\"\n",
    "OUTPUT_DIR = r\"E:\\texas_precip\"\n",
    "\n",
    "# (filename_code, duration_minutes) — matches NOAA naming convention\n",
    "DURATIONS = [\n",
    "    (\"05m\",  5),\n",
    "    (\"10m\",  10),\n",
    "    (\"15m\",  15),\n",
    "    (\"30m\",  30),\n",
    "    (\"45m\",  45),\n",
    "    (\"60m\",  60),\n",
    "    (\"02h\",  120),\n",
    "    (\"03h\",  180),\n",
    "    (\"06h\",  360),\n",
    "    (\"12h\",  720),\n",
    "    (\"24h\",  1440),\n",
    "]\n",
    "\n",
    "os.makedirs(OUTPUT_DIR, exist_ok=True)\n",
    "\n",
    "\n",
    "def download_and_unzip(url, dest_dir):\n",
    "    filename = url.split(\"/\")[-1]\n",
    "    filepath = os.path.join(dest_dir, filename)\n",
    "\n",
    "    print(f\"Downloading {filename} ...\", end=\" \", flush=True)\n",
    "    try:\n",
    "        urllib.request.urlretrieve(url, filepath)\n",
    "        print(\"OK\", end=\" \", flush=True)\n",
    "    except Exception as e:\n",
    "        print(f\"FAILED ({e})\")\n",
    "        return\n",
    "\n",
    "    try:\n",
    "        with zipfile.ZipFile(filepath, \"r\") as z:\n",
    "            z.extractall(dest_dir)\n",
    "        os.remove(filepath)\n",
    "        print(\"-> unzipped\")\n",
    "    except zipfile.BadZipFile:\n",
    "        print(\"-> not a valid zip, kept as-is\")\n",
    "\n",
    "\n",
    "for code, dur_min in DURATIONS:\n",
    "    filename = f\"tx1yr{code}a.zip\"\n",
    "    url = BASE_URL + filename\n",
    "    download_and_unzip(url, OUTPUT_DIR)\n",
    "\n",
    "print(f\"\\nDone. Files saved to: {os.path.abspath(OUTPUT_DIR)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b01f5b96",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Convert NOAA Atlas 14 Texas 1-year ASC precipitation files\n",
    "(all durations up to 24 hours) to a single NetCDF4 file.\n",
    "\n",
    "- Input : tx1yr{dur}a.asc files in INPUT_DIR\n",
    "- Output: tx_1yr_intensity.nc in OUTPUT_DIR\n",
    "- Structure: single variable intensity(duration_minutes, lat, lon)\n",
    "- Values: converted from integer (thousandths of inches depth)\n",
    "          to intensity in inches/hour  (depth_in / duration_hr)\n",
    "- Missing: nodata -> NaN\n",
    "\"\"\"\n",
    "\n",
    "import os\n",
    "import re\n",
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "\n",
    "# --- Configuration ---\n",
    "INPUT_DIR = r\"E:\\texas_precip\"\n",
    "OUTPUT_NC = r\"E:\\texas_precip\\tx_1yr_intensity.nc\"\n",
    "\n",
    "# (filename_code, duration_minutes)\n",
    "DURATIONS = [\n",
    "    (\"05m\",    5),\n",
    "    (\"10m\",   10),\n",
    "    (\"15m\",   15),\n",
    "    (\"30m\",   30),\n",
    "    (\"45m\",   45),\n",
    "    (\"60m\",   60),\n",
    "    (\"02h\",  120),\n",
    "    (\"03h\",  180),\n",
    "    (\"06h\",  360),\n",
    "    (\"12h\",  720),\n",
    "    (\"24h\", 1440),\n",
    "]\n",
    "\n",
    "NODATA_RAW = -9\n",
    "SCALE      = 1000.0   # values stored as thousandths of inches\n",
    "\n",
    "\n",
    "def parse_asc_header(f):\n",
    "    header = {}\n",
    "    for _ in range(6):\n",
    "        line = f.readline().strip().split()\n",
    "        header[line[0].lower()] = float(line[1])\n",
    "    return header\n",
    "\n",
    "\n",
    "def read_asc(filepath):\n",
    "    with open(filepath, \"r\") as f:\n",
    "        hdr = parse_asc_header(f)\n",
    "        data = np.loadtxt(f, dtype=np.float32)\n",
    "    return hdr, data\n",
    "\n",
    "\n",
    "def to_intensity(raw, nodata, dur_min):\n",
    "    \"\"\"Convert raw integer grid to intensity in inches/hour.\"\"\"\n",
    "    data = raw.astype(np.float32)\n",
    "    data[data == nodata]     = np.nan\n",
    "    data[data == NODATA_RAW] = np.nan\n",
    "    depth_in = data / SCALE\n",
    "    dur_hr   = dur_min / 60.0\n",
    "    return depth_in / dur_hr\n",
    "\n",
    "\n",
    "# --- Discover files ---\n",
    "found = {}\n",
    "for code, dur_min in DURATIONS:\n",
    "    fname = f\"tx1yr{code}a.asc\"\n",
    "    fpath = os.path.join(INPUT_DIR, fname)\n",
    "    if os.path.exists(fpath):\n",
    "        found[dur_min] = (code, fpath)\n",
    "    else:\n",
    "        print(f\"WARNING: {fname} not found in {INPUT_DIR}\")\n",
    "\n",
    "durs_to_use = sorted(found.keys())\n",
    "if not durs_to_use:\n",
    "    raise FileNotFoundError(f\"No matching ASC files found in {INPUT_DIR}\")\n",
    "\n",
    "print(f\"Found {len(durs_to_use)} ASC files for durations (min): {durs_to_use}\")\n",
    "\n",
    "# --- Read first file to get grid geometry ---\n",
    "_, first_path = found[durs_to_use[0]]\n",
    "hdr0, _ = read_asc(first_path)\n",
    "ncols  = int(hdr0[\"ncols\"])\n",
    "nrows  = int(hdr0[\"nrows\"])\n",
    "xll    = hdr0[\"xllcorner\"]\n",
    "yll    = hdr0[\"yllcorner\"]\n",
    "cell   = hdr0[\"cellsize\"]\n",
    "nodata = hdr0.get(\"nodata_value\", NODATA_RAW)\n",
    "\n",
    "lons = xll + (np.arange(ncols) + 0.5) * cell\n",
    "lats = yll + (np.arange(nrows - 1, -1, -1) + 0.5) * cell\n",
    "\n",
    "# --- Write NetCDF ---\n",
    "print(f\"Writing {OUTPUT_NC} ...\")\n",
    "ds = nc.Dataset(OUTPUT_NC, \"w\", format=\"NETCDF4\")\n",
    "\n",
    "ds.title        = \"NOAA Atlas 14 Texas 1-Year Precipitation Intensity by Duration\"\n",
    "ds.institution  = \"NOAA/NWS Hydrometeorological Design Studies Center\"\n",
    "ds.source       = \"NOAA Atlas 14 Volume 11\"\n",
    "ds.Conventions  = \"CF-1.8\"\n",
    "ds.return_period = \"1 year\"\n",
    "ds.crs          = \"NAD83 geographic (EPSG:4269)\"\n",
    "\n",
    "# Dimensions\n",
    "ds.createDimension(\"duration_minutes\", len(durs_to_use))\n",
    "ds.createDimension(\"lat\", nrows)\n",
    "ds.createDimension(\"lon\", ncols)\n",
    "\n",
    "# Coordinate variables\n",
    "dur_var = ds.createVariable(\"duration_minutes\", \"i4\", (\"duration_minutes\",))\n",
    "dur_var.units     = \"minutes\"\n",
    "dur_var.long_name = \"storm duration\"\n",
    "dur_var[:]        = np.array(durs_to_use, dtype=np.int32)\n",
    "\n",
    "lat_var = ds.createVariable(\"lat\", \"f4\", (\"lat\",))\n",
    "lat_var.units         = \"degrees_north\"\n",
    "lat_var.standard_name = \"latitude\"\n",
    "lat_var[:]            = lats\n",
    "\n",
    "lon_var = ds.createVariable(\"lon\", \"f4\", (\"lon\",))\n",
    "lon_var.units         = \"degrees_east\"\n",
    "lon_var.standard_name = \"longitude\"\n",
    "lon_var[:]            = lons\n",
    "\n",
    "# Single intensity variable\n",
    "intensity_var = ds.createVariable(\n",
    "    \"intensity\", \"f4\", (\"duration_minutes\", \"lat\", \"lon\"),\n",
    "    fill_value=np.nan, zlib=True, complevel=4\n",
    ")\n",
    "intensity_var.long_name     = \"1-year precipitation intensity\"\n",
    "intensity_var.units         = \"inches/hour\"\n",
    "intensity_var.return_period = \"1 year\"\n",
    "intensity_var.grid_mapping  = \"crs\"\n",
    "intensity_var.coordinates   = \"duration_minutes lat lon\"\n",
    "\n",
    "# Fill intensity array\n",
    "for i, dur_min in enumerate(durs_to_use):\n",
    "    code, fpath = found[dur_min]\n",
    "    dur_hr = dur_min / 60.0\n",
    "    print(f\"  Processing {code} ({dur_min} min) ...\", end=\" \", flush=True)\n",
    "    hdr, raw = read_asc(fpath)\n",
    "    intensity_var[i, :, :] = to_intensity(raw, nodata, dur_min)\n",
    "    print(\"OK\")\n",
    "\n",
    "# CRS variable\n",
    "crs_var = ds.createVariable(\"crs\", \"i4\")\n",
    "crs_var.grid_mapping_name           = \"latitude_longitude\"\n",
    "crs_var.longitude_of_prime_meridian = 0.0\n",
    "crs_var.semi_major_axis             = 6378137.0\n",
    "crs_var.inverse_flattening          = 298.257222\n",
    "crs_var.crs_wkt = (\n",
    "    'GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",'\n",
    "    'SPHEROID[\"GRS 1980\",6378137,298.257222101]],'\n",
    "    'PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]'\n",
    ")\n",
    "\n",
    "ds.close()\n",
    "print(f\"\\nDone. Output: {OUTPUT_NC}\")\n",
    "print(f\"  Grid      : {nrows} rows x {ncols} cols\")\n",
    "print(f\"  Durations : {durs_to_use} minutes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dfabb12",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Combine tx_1yr_intensity.nc and tx_5min_intensity.nc into a single NetCDF.\n",
    "\n",
    "Output variables:\n",
    "  intensity_by_duration(duration_minutes, lat, lon)\n",
    "      - 1-year return period, durations 5-1440 min, units: inches/hour\n",
    "\n",
    "  intensity_by_return_period(return_period, lat, lon)\n",
    "      - 5-minute duration, return periods 1-1000 yr, units: inches/hour\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "\n",
    "INPUT_1YR  = r\"E:\\texas_precip\\tx_1yr_intensity.nc\"\n",
    "INPUT_5MIN = r\"E:\\texas_precip\\tx_5min_intensity.nc\"\n",
    "OUTPUT_NC  = r\"E:\\texas_precip\\tx_precip_intensity.nc\"\n",
    "\n",
    "print(\"Opening source files...\")\n",
    "src1 = nc.Dataset(INPUT_1YR,  \"r\")\n",
    "src5 = nc.Dataset(INPUT_5MIN, \"r\")\n",
    "\n",
    "# --- Verify grids match ---\n",
    "assert np.allclose(src1[\"lat\"][:], src5[\"lat\"][:]), \"lat grids do not match!\"\n",
    "assert np.allclose(src1[\"lon\"][:], src5[\"lon\"][:]), \"lon grids do not match!\"\n",
    "lats = src1[\"lat\"][:]\n",
    "lons = src1[\"lon\"][:]\n",
    "nrows, ncols = len(lats), len(lons)\n",
    "print(f\"Grid: {nrows} rows x {ncols} cols — grids match OK\")\n",
    "\n",
    "# --- Create output ---\n",
    "print(f\"Writing {OUTPUT_NC} ...\")\n",
    "ds = nc.Dataset(OUTPUT_NC, \"w\", format=\"NETCDF4\")\n",
    "\n",
    "ds.title       = \"NOAA Atlas 14 Texas Precipitation Intensity\"\n",
    "ds.institution = \"NOAA/NWS Hydrometeorological Design Studies Center\"\n",
    "ds.source      = \"NOAA Atlas 14 Volume 11\"\n",
    "ds.Conventions = \"CF-1.8\"\n",
    "ds.crs         = \"NAD83 geographic (EPSG:4269)\"\n",
    "ds.contents    = (\n",
    "    \"intensity_by_duration: 1-year return period, variable duration; \"\n",
    "    \"intensity_by_return_period: 5-minute duration, variable return period\"\n",
    ")\n",
    "\n",
    "# Dimensions\n",
    "dur_vals = src1[\"duration_minutes\"][:]\n",
    "rp_vals  = src5[\"return_period\"][:]\n",
    "\n",
    "ds.createDimension(\"duration_minutes\", len(dur_vals))\n",
    "ds.createDimension(\"return_period\",    len(rp_vals))\n",
    "ds.createDimension(\"lat\",  nrows)\n",
    "ds.createDimension(\"lon\",  ncols)\n",
    "\n",
    "# Coordinate variables\n",
    "dur_var = ds.createVariable(\"duration_minutes\", \"i4\", (\"duration_minutes\",))\n",
    "dur_var.units     = \"minutes\"\n",
    "dur_var.long_name = \"storm duration\"\n",
    "dur_var[:]        = dur_vals\n",
    "\n",
    "rp_var = ds.createVariable(\"return_period\", \"i4\", (\"return_period\",))\n",
    "rp_var.units     = \"years\"\n",
    "rp_var.long_name = \"return period\"\n",
    "rp_var[:]        = rp_vals\n",
    "\n",
    "lat_var = ds.createVariable(\"lat\", \"f4\", (\"lat\",))\n",
    "lat_var.units         = \"degrees_north\"\n",
    "lat_var.standard_name = \"latitude\"\n",
    "lat_var[:]            = lats\n",
    "\n",
    "lon_var = ds.createVariable(\"lon\", \"f4\", (\"lon\",))\n",
    "lon_var.units         = \"degrees_east\"\n",
    "lon_var.standard_name = \"longitude\"\n",
    "lon_var[:]            = lons\n",
    "\n",
    "# Variable 1: 1-year intensity by duration\n",
    "print(\"  Copying intensity_by_duration ...\")\n",
    "v1 = ds.createVariable(\n",
    "    \"intensity_by_duration\", \"f4\", (\"duration_minutes\", \"lat\", \"lon\"),\n",
    "    fill_value=np.nan, zlib=True, complevel=4\n",
    ")\n",
    "v1.long_name     = \"1-year precipitation intensity by storm duration\"\n",
    "v1.units         = \"inches/hour\"\n",
    "v1.return_period = \"1 year\"\n",
    "v1.grid_mapping  = \"crs\"\n",
    "v1[:]            = src1[\"intensity\"][:]\n",
    "\n",
    "# Variable 2: 5-minute intensity by return period\n",
    "print(\"  Copying intensity_by_return_period ...\")\n",
    "v2 = ds.createVariable(\n",
    "    \"intensity_by_return_period\", \"f4\", (\"return_period\", \"lat\", \"lon\"),\n",
    "    fill_value=np.nan, zlib=True, complevel=4\n",
    ")\n",
    "v2.long_name      = \"5-minute precipitation intensity by return period\"\n",
    "v2.units          = \"inches/hour\"\n",
    "v2.duration_minutes = 5\n",
    "v2.grid_mapping   = \"crs\"\n",
    "v2[:]             = src5[\"intensity\"][:]\n",
    "\n",
    "# CRS variable\n",
    "crs_var = ds.createVariable(\"crs\", \"i4\")\n",
    "crs_var.grid_mapping_name           = \"latitude_longitude\"\n",
    "crs_var.longitude_of_prime_meridian = 0.0\n",
    "crs_var.semi_major_axis             = 6378137.0\n",
    "crs_var.inverse_flattening          = 298.257222\n",
    "crs_var.crs_wkt = (\n",
    "    'GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",'\n",
    "    'SPHEROID[\"GRS 1980\",6378137,298.257222101]],'\n",
    "    'PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]'\n",
    ")\n",
    "\n",
    "src1.close()\n",
    "src5.close()\n",
    "ds.close()\n",
    "\n",
    "print(f\"\\nDone. Output: {OUTPUT_NC}\")\n",
    "print(f\"  intensity_by_duration      : {len(dur_vals)} durations (min): {list(dur_vals)}\")\n",
    "print(f\"  intensity_by_return_period : {len(rp_vals)} return periods (yr): {list(rp_vals)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8df66c4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Combine tx_1yr_intensity.nc and tx_5min_intensity.nc into a single NetCDF.\n",
    "\n",
    "Output variables:\n",
    "  intensity_by_duration(duration_minutes, lat, lon)\n",
    "      - 1-year return period, durations 5-1440 min, units: inches/hour\n",
    "\n",
    "  intensity_by_return_period(return_period, lat, lon)\n",
    "      - 5-minute duration, return periods 1-1000 yr, units: inches/hour\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "\n",
    "INPUT_1YR  = r\"E:\\texas_precip\\tx_1yr_intensity.nc\"\n",
    "INPUT_5MIN = r\"E:\\texas_precip\\tx_5min_intensity.nc\"\n",
    "OUTPUT_NC  = r\"E:\\texas_precip\\tx_precip_intensity.nc\"\n",
    "\n",
    "print(\"Opening source files...\")\n",
    "src1 = nc.Dataset(INPUT_1YR,  \"r\")\n",
    "src5 = nc.Dataset(INPUT_5MIN, \"r\")\n",
    "\n",
    "# --- Verify grids match ---\n",
    "assert np.allclose(src1[\"lat\"][:], src5[\"lat\"][:]), \"lat grids do not match!\"\n",
    "assert np.allclose(src1[\"lon\"][:], src5[\"lon\"][:]), \"lon grids do not match!\"\n",
    "lats = src1[\"lat\"][:]\n",
    "lons = src1[\"lon\"][:]\n",
    "nrows, ncols = len(lats), len(lons)\n",
    "print(f\"Grid: {nrows} rows x {ncols} cols — grids match OK\")\n",
    "\n",
    "# --- Create output ---\n",
    "print(f\"Writing {OUTPUT_NC} ...\")\n",
    "ds = nc.Dataset(OUTPUT_NC, \"w\", format=\"NETCDF4\")\n",
    "\n",
    "ds.title       = \"NOAA Atlas 14 Texas Precipitation Intensity\"\n",
    "ds.institution = \"NOAA/NWS Hydrometeorological Design Studies Center\"\n",
    "ds.source      = \"NOAA Atlas 14 Volume 11\"\n",
    "ds.Conventions = \"CF-1.8\"\n",
    "ds.crs         = \"WGS84 geographic (EPSG:4326)\"\n",
    "ds.contents    = (\n",
    "    \"intensity_by_duration: 1-year return period, variable duration; \"\n",
    "    \"intensity_by_return_period: 5-minute duration, variable return period\"\n",
    ")\n",
    "\n",
    "# Dimensions\n",
    "dur_vals = src1[\"duration_minutes\"][:]\n",
    "rp_vals  = src5[\"return_period\"][:]\n",
    "\n",
    "ds.createDimension(\"duration_minutes\", len(dur_vals))\n",
    "ds.createDimension(\"return_period\",    len(rp_vals))\n",
    "ds.createDimension(\"lat\",  nrows)\n",
    "ds.createDimension(\"lon\",  ncols)\n",
    "\n",
    "# Coordinate variables\n",
    "dur_var = ds.createVariable(\"duration_minutes\", \"i4\", (\"duration_minutes\",))\n",
    "dur_var.units     = \"minutes\"\n",
    "dur_var.long_name = \"storm duration\"\n",
    "dur_var[:]        = dur_vals\n",
    "\n",
    "rp_var = ds.createVariable(\"return_period\", \"i4\", (\"return_period\",))\n",
    "rp_var.units     = \"years\"\n",
    "rp_var.long_name = \"return period\"\n",
    "rp_var[:]        = rp_vals\n",
    "\n",
    "lat_var = ds.createVariable(\"lat\", \"f4\", (\"lat\",))\n",
    "lat_var.units         = \"degrees_north\"\n",
    "lat_var.standard_name = \"latitude\"\n",
    "lat_var[:]            = lats\n",
    "\n",
    "lon_var = ds.createVariable(\"lon\", \"f4\", (\"lon\",))\n",
    "lon_var.units         = \"degrees_east\"\n",
    "lon_var.standard_name = \"longitude\"\n",
    "lon_var[:]            = lons\n",
    "\n",
    "# Variable 1: 1-year intensity by duration\n",
    "print(\"  Copying intensity_by_duration ...\")\n",
    "v1 = ds.createVariable(\n",
    "    \"intensity_by_duration\", \"f4\", (\"duration_minutes\", \"lat\", \"lon\"),\n",
    "    fill_value=np.nan, zlib=True, complevel=4\n",
    ")\n",
    "v1.long_name     = \"1-year precipitation intensity by storm duration\"\n",
    "v1.units         = \"inches/hour\"\n",
    "v1.return_period = \"1 year\"\n",
    "v1.grid_mapping  = \"crs\"\n",
    "v1[:]            = src1[\"intensity\"][:]\n",
    "\n",
    "# Variable 2: 5-minute intensity by return period\n",
    "print(\"  Copying intensity_by_return_period ...\")\n",
    "v2 = ds.createVariable(\n",
    "    \"intensity_by_return_period\", \"f4\", (\"return_period\", \"lat\", \"lon\"),\n",
    "    fill_value=np.nan, zlib=True, complevel=4\n",
    ")\n",
    "v2.long_name      = \"5-minute precipitation intensity by return period\"\n",
    "v2.units          = \"inches/hour\"\n",
    "v2.duration_minutes = 5\n",
    "v2.grid_mapping   = \"crs\"\n",
    "v2[:]             = src5[\"intensity\"][:]\n",
    "\n",
    "# CRS variable\n",
    "crs_var = ds.createVariable(\"crs\", \"i4\")\n",
    "crs_var.grid_mapping_name           = \"latitude_longitude\"\n",
    "crs_var.longitude_of_prime_meridian = 0.0\n",
    "crs_var.semi_major_axis             = 6378137.0\n",
    "crs_var.inverse_flattening          = 298.257223563\n",
    "crs_var.epsg_code                   = \"EPSG:4326\"\n",
    "crs_var.crs_wkt = (\n",
    "    'GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",'\n",
    "    'SPHEROID[\"WGS 84\",6378137,298.257223563]],'\n",
    "    'PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],'\n",
    "    'AUTHORITY[\"EPSG\",\"4326\"]]'\n",
    ")\n",
    "\n",
    "src1.close()\n",
    "src5.close()\n",
    "ds.close()\n",
    "\n",
    "print(f\"\\nDone. Output: {OUTPUT_NC}\")\n",
    "print(f\"  intensity_by_duration      : {len(dur_vals)} durations (min): {list(dur_vals)}\")\n",
    "print(f\"  intensity_by_return_period : {len(rp_vals)} return periods (yr): {list(rp_vals)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71bb1f2e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}