{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "b52348f6", "metadata": {}, "outputs": [], "source": [ "# USES mrms-grib env" ] }, { "cell_type": "code", "execution_count": 57, "id": "0547d8cf", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: total: 1.73 s\n", "Wall time: 4.67 s\n" ] } ], "source": [ "%%time\n", "# --- Get the list of intensities from Precipitation netCDF\n", "# --- given a field name and value from a polygon dataset\n", "# --- assumes that the dataset is EPSG:5070\n", "# --- and netCDF is EPSG:4326\n", "\n", "import netCDF4 as nc\n", "import numpy as np\n", "from shapely.geometry import shape, Point\n", "from shapely.prepared import prep\n", "import geopandas as gpd\n", "\n", "import warnings\n", "import numpy as np\n", "\n", "\n", "nc_path=r'E:\\texas_precip\\tx_precip_intensity.nc'\n", "\n", "fgb_path = 'https://web.corral.tacc.utexas.edu/ras2fim/huc12/huc12_texas_ar_5070.fgb'\n", "str_field_to_search = 'huc12'\n", "str_id_to_find = '120301020702'\n", "\n", "#fgb_path = 'https://web.corral.tacc.utexas.edu/ras2fim/nextgen/divides_nextgen_12.fgb'\n", "#str_field_to_search = 'divide_id'\n", "#str_id_to_find = 'cat-2402391'\n", "\n", "# ----------------------------\n", "def fn_get_peak_intensity_pixel(polygon, nc_path,\n", " return_period_idx=None, agg='max'):\n", " \"\"\"\n", " Find the pixel with the highest intensity_by_return_period value\n", " that intersects a given polygon.\n", "\n", " Parameters\n", " ----------\n", " polygon : shapely.geometry.Polygon or dict\n", " A Shapely polygon or GeoJSON-style geometry dict defining the area of interest.\n", " nc_path : str\n", " Path to the NetCDF file.\n", " return_period_idx : int or None\n", " Index into the return_period dimension (0–9) to use.\n", " If None, uses the maximum across ALL return periods.\n", " agg : str\n", " How to aggregate across return periods when return_period_idx is None.\n", " 'max' = peak across all return periods (default).\n", "\n", " Returns\n", " -------\n", " dict with keys:\n", " lat : float - latitude of peak pixel\n", " lon : float - longitude of peak pixel\n", " lat_idx : int - row index in the grid\n", " lon_idx : int - col index in the grid\n", " intensity : float - intensity value (inches/hour)\n", " return_period : int - return period (years) at which peak occurs\n", " return_period_idx : int - index of that return period\n", " n_pixels_in_polygon : int - total pixels found inside polygon\n", " \"\"\"\n", " # Accept GeoJSON dict as well as Shapely geometry\n", " if isinstance(polygon, dict):\n", " polygon = shape(polygon)\n", "\n", " with nc.Dataset(nc_path, 'r') as ds:\n", " lats = ds.variables['lat'][:] # (1282,)\n", " lons = ds.variables['lon'][:] # (1579,)\n", " return_periods = ds.variables['return_period'][:] # (10,)\n", "\n", " # --- 1. Bounding-box pre-filter (fast) ---\n", " minx, miny, maxx, maxy = polygon.bounds\n", "\n", " lat_mask = (lats >= miny) & (lats <= maxy)\n", " lon_mask = (lons >= minx) & (lons <= maxx)\n", "\n", " lat_indices = np.where(lat_mask)[0]\n", " lon_indices = np.where(lon_mask)[0]\n", "\n", " if lat_indices.size == 0 or lon_indices.size == 0:\n", " raise ValueError(\"Polygon does not overlap the NetCDF grid extent.\")\n", "\n", " lat_slice = slice(lat_indices[0], lat_indices[-1] + 1)\n", " lon_slice = slice(lon_indices[0], lon_indices[-1] + 1)\n", "\n", " # Load only the bounding-box subset of the intensity grid\n", " # shape: (10, n_lat_bb, n_lon_bb)\n", " if return_period_idx is not None:\n", " # Single return period — load one slice\n", " intensity_bb = ds.variables['intensity_by_return_period'][\n", " return_period_idx, lat_slice, lon_slice\n", " ] # (n_lat_bb, n_lon_bb)\n", " intensity_bb = intensity_bb[np.newaxis, ...] # keep 3-D for uniform code\n", " rp_indices_loaded = [return_period_idx]\n", " else:\n", " intensity_bb = ds.variables['intensity_by_return_period'][\n", " :, lat_slice, lon_slice\n", " ] # (10, n_lat_bb, n_lon_bb)\n", " rp_indices_loaded = list(range(len(return_periods)))\n", "\n", " lats_bb = lats[lat_slice] # (n_lat_bb,)\n", " lons_bb = lons[lon_slice] # (n_lon_bb,)\n", "\n", " # --- 2. Point-in-polygon test (prepared geometry for speed) ---\n", " prepared_poly = prep(polygon)\n", "\n", " # Build coordinate grids for the bounding box\n", " lon_grid, lat_grid = np.meshgrid(lons_bb, lats_bb) # both (n_lat_bb, n_lon_bb)\n", "\n", " n_lat_bb, n_lon_bb = lat_grid.shape\n", " inside = np.zeros((n_lat_bb, n_lon_bb), dtype=bool)\n", "\n", " for i in range(n_lat_bb):\n", " for j in range(n_lon_bb):\n", " inside[i, j] = prepared_poly.contains(Point(lon_grid[i, j], lat_grid[i, j]))\n", "\n", " n_pixels = int(inside.sum())\n", " if n_pixels == 0:\n", " raise ValueError(\"No grid pixels found inside the polygon. \"\n", " \"The polygon may be smaller than one grid cell — \"\n", " \"consider using intersects() instead of contains().\")\n", "\n", " # --- 3. Find peak intensity across pixels inside the polygon ---\n", " # Mask pixels outside polygon\n", " inside_3d = np.broadcast_to(inside[np.newaxis, :, :], intensity_bb.shape)\n", " intensity_masked = np.where(inside_3d, intensity_bb, np.nan)\n", "\n", " # Max across all loaded return periods → best pixel\n", " intensity_max_over_rp = np.nanmax(intensity_masked, axis=0) # (n_lat_bb, n_lon_bb)\n", "\n", " # Pixel with highest intensity\n", " flat_idx = np.nanargmax(intensity_max_over_rp)\n", " best_i, best_j = np.unravel_index(flat_idx, (n_lat_bb, n_lon_bb))\n", "\n", " # Which return period gave that max?\n", " best_rp_local = int(np.nanargmax(intensity_masked[:, best_i, best_j]))\n", " best_rp_idx = rp_indices_loaded[best_rp_local]\n", "\n", " # Map back to global grid indices\n", " global_lat_idx = int(lat_indices[0] + best_i)\n", " global_lon_idx = int(lon_indices[0] + best_j)\n", "\n", " return {\n", " 'lat': float(lats_bb[best_i]),\n", " 'lon': float(lons_bb[best_j]),\n", " 'lat_idx': global_lat_idx,\n", " 'lon_idx': global_lon_idx,\n", " 'intensity': float(intensity_max_over_rp[best_i, best_j]),\n", " 'return_period': int(return_periods[best_rp_idx]),\n", " 'return_period_idx': best_rp_idx,\n", " 'n_pixels_in_polygon': n_pixels,\n", " }\n", "# ----------------------------\n", "\n", "# --------------------------\n", "def fn_get_polygon_by_value(str_fieldname: str, str_value: str, fgb_path):\n", " \"\"\"\n", " Given a field name and value, return the corresponding Shapely polygon.\n", "\n", " Parameters\n", " ----------\n", " str_fieldname : str\n", " The field/column name to filter on.\n", " str_value : str\n", " The value to match.\n", " fgb_path : str\n", " Path to the local FlatGeobuf file.\n", "\n", " Returns\n", " -------\n", " shapely.geometry.Polygon or MultiPolygon\n", " \"\"\"\n", " import fiona\n", " from shapely.geometry import shape\n", "\n", " with fiona.open(fgb_path) as src:\n", "\n", " # validate field name\n", " valid_fields = list(src.schema['properties'].keys())\n", " if str_fieldname not in valid_fields:\n", " raise KeyError(\n", " f\"Field '{str_fieldname}' not found. \"\n", " f\"Available fields: {valid_fields}\"\n", " )\n", "\n", " for feature in src:\n", " if str(feature['properties'][str_fieldname]) == str(str_value):\n", " return shape(feature['geometry'])\n", "\n", " return None\n", "# --------------------------\n", "\n", "# ----------------------\n", "def fn_get_intensity_by_duration(lat_idx: int, lon_idx: int, nc_path):\n", " with nc.Dataset(nc_path, 'r') as ds:\n", " intensities = ds.variables['intensity_by_duration'][:, lat_idx, lon_idx].tolist()\n", " return sorted(intensities)\n", "# ----------------------\n", "\n", "# ----------------------\n", "def fn_get_intensity_by_return_period(lat_idx: int, lon_idx: int, nc_path):\n", " with nc.Dataset(nc_path, 'r') as ds:\n", " intensities = ds.variables['intensity_by_return_period'][:, lat_idx, lon_idx].tolist()\n", " return sorted(intensities)\n", "# ----------------------\n", "\n", "\n", "shape_poly_5070 = fn_get_polygon_by_value(str_field_to_search, str_id_to_find, fgb_path)\n", "\n", "if shape_poly_5070 is None:\n", " print(f\"Field:{str_field_to_search} Value:{str_id_to_find} No feature found.\")\n", "else:\n", " shape_poly_4326 = (gpd.GeoDataFrame(geometry=[shape_poly_5070], crs='EPSG:5070').to_crs('EPSG:4326').geometry.iloc[0])\n", " with warnings.catch_warnings():\n", " warnings.simplefilter(\"ignore\", RuntimeWarning)\n", " dict_peak_pixel_result = fn_get_peak_intensity_pixel(shape_poly_4326, nc_path)\n", "\n", " list_intensity_in_hr_1year = fn_get_intensity_by_duration(\n", " dict_peak_pixel_result['lat_idx'],\n", " dict_peak_pixel_result['lon_idx'],\n", " nc_path)\n", "\n", " list_intensity_by_return_period_in_hr_5min_duration = fn_get_intensity_by_return_period(\n", " dict_peak_pixel_result['lat_idx'],\n", " dict_peak_pixel_result['lon_idx'],\n", " nc_path)\n", "\n", " list_in_hr_intensities = sorted(set(list_intensity_in_hr_1year + list_intensity_by_return_period_in_hr_5min_duration))\n", "\n", " list_mm_hr_intensities = [round(x * 25.4, 1) for x in list_in_hr_intensities]" ] }, { "cell_type": "code", "execution_count": 58, "id": "6ec3106f-8e07-45c4-b9cb-701fa5087c74", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[3.4,\n", " 5.9,\n", " 10.0,\n", " 17.0,\n", " 23.1,\n", " 37.9,\n", " 58.3,\n", " 83.8,\n", " 100.7,\n", " 125.9,\n", " 147.5,\n", " 182.6,\n", " 211.8,\n", " 252.1,\n", " 282.9,\n", " 314.2,\n", " 347.5,\n", " 392.6,\n", " 427.9]" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_mm_hr_intensities" ] }, { "cell_type": "code", "execution_count": 59, "id": "a7aef075-9f44-4301-94d5-d4fc8d3ebbc2", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gpd.GeoDataFrame(geometry=[shape_poly_4326], crs='EPSG:4326').explore()" ] }, { "cell_type": "code", "execution_count": null, "id": "6791fabe-633d-4181-a43a-ac51b981b332", "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.11.14" } }, "nbformat": 4, "nbformat_minor": 5 }