diff --git a/requirements.dev.txt b/requirements.dev.txt deleted file mode 100644 index 486c66c..0000000 --- a/requirements.dev.txt +++ /dev/null @@ -1,19 +0,0 @@ -Cython==3.0.10 -dask==2024.4.1 -gcsfs==2024.3.1 -ipykernel==6.29.4 -Jinja2==3.1.3 -matplotlib==3.8.4 -MetPy==1.6.2 -numpy==1.26.4 -xarray==2024.3.0 -zarr==2.17.2 -# - testing -coverage==7.4.4 -pytest==7.4.2 -pytest-cov==4.1.0 -# - linting -ruff==0.3.7 -black==24.4.0 -isort==5.13.2 -jupyter-black==0.3.4 \ No newline at end of file diff --git a/scratchpad.ipynb b/scratchpad.ipynb deleted file mode 100644 index ea3bda7..0000000 --- a/scratchpad.ipynb +++ /dev/null @@ -1,299 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Ignoring index file 'data/hrrr.t00z.wrfprsf00.grib2.9093e.idx' incompatible with GRIB file\n" - ] - } - ], - "source": [ - "from __future__ import annotations\n", - "\n", - "\n", - "import datetime\n", - "from typing import Any\n", - "\n", - "import numpy as np\n", - "import xarray as xr\n", - "import metpy.calc as mpcalc\n", - "from metpy.units import units\n", - "\n", - "import src.nzthermo as nzt \n", - "# from src.nzthermo._core import nanshift\n", - "\n", - "\n", - "\n", - "Pa = units.pascal\n", - "K = units.kelvin\n", - "\n", - "isobaric = xr.open_dataset(\n", - " \"data/hrrr.t00z.wrfprsf00.grib2\",\n", - " engine=\"cfgrib\",\n", - " backend_kwargs={\"filter_by_keys\": {\"typeOfLevel\": \"isobaricInhPa\"}},\n", - ")\n", - "\n", - "\n", - "def get_data(**sel: Any):\n", - " ds = isobaric.sel(**sel)\n", - " T = ds[\"t\"].to_numpy() # (K) (Z, Y, X)\n", - " Z, Y, X = T.shape\n", - " N = Y * X\n", - " T = T.reshape(Z, N).transpose() # (N, Z)\n", - " P = ds[\"isobaricInhPa\"].to_numpy().astype(np.float32) * 100.0 # (Pa)\n", - " Q = ds[\"q\"].to_numpy() # (kg/kg) (Z, Y, X)\n", - " Q = Q.reshape(Z, N).transpose() # (N, Z)\n", - " Td = nzt.dewpoint_from_specific_humidity(P, Q)\n", - " lat = ds[\"latitude\"].to_numpy()\n", - " lon = (ds[\"longitude\"].to_numpy() + 180) % 360 - 180\n", - " timestamp = datetime.datetime.fromisoformat(ds[\"time\"].to_numpy().astype(str).item())\n", - " extent = [lon.min(), lon.max(), lat.min(), lat.max()]\n", - "\n", - " return (P, T, Td), (Z, Y, X), (lat, lon, timestamp, extent)\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "mpcalc.most_unstable_parcel # \n", - "(P, T, Td), (Z, Y, X), (latitude, longitude, timestamp, extent) = get_data(\n", - " x = slice(None, None, 45), \n", - " y = slice(None, None, 45),\n", - ")\n", - "\n", - "\n", - "from src.nzthermo import equivalent_potential_temperature, parcel_profile_with_lcl\n", - "\n", - "\n", - "\n", - "\n", - "def most_unstable_parcel_index(\n", - " pressure,\n", - " temperature,\n", - " dewpoint,\n", - " /,\n", - " depth: float = 30_000.0,\n", - " height: float | None = None,\n", - " bottom: float | None = None,\n", - "):\n", - " if height is not None:\n", - " raise NotImplementedError(\"height argument is not implemented\")\n", - " \n", - " pressure = np.atleast_2d(pressure)\n", - " p0 = (pressure[:, 0] if bottom is None else np.asarray(bottom))#.reshape(-1, 1)\n", - " top = p0 - depth\n", - " mask, = np.nonzero(nzt._ufunc.between_or_close(pressure, top, p0).astype(np.bool_).squeeze())\n", - " t_layer = temperature[:, mask]\n", - " td_layer = dewpoint[:, mask]\n", - " p_layer = pressure[:, mask]\n", - " theta_e = equivalent_potential_temperature(p_layer, t_layer, td_layer)\n", - " return np.argmax(theta_e, axis=1)\n", - "\n", - "\n", - "np.testing.assert_array_equal(\n", - " most_unstable_parcel_index(P, T, Td), \n", - " [mpcalc.most_unstable_parcel(P*Pa, T[i]*K, Td[i]*K)[-1] for i in range(Y*X)]\n", - ")\n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# idx = most_unstable_parcel(P, T, Td)\n", - "\n", - "# mask = ~(idx[:, None] >= np.arange(Z)[None, :])\n", - "\n", - "# # mask.argmin(1) <= np.arange(Z)\n", - "\n", - "# P_new = np.full_like(mask, np.nan)\n", - "\n", - "\n", - "\n", - "\n", - "# # # np.argmin(~mask, axis=1)\n", - "# # # P_new\n", - "\n", - "# # np.nonzero(mask)\n", - "\n", - "# np.where(mask, P, np.nan)\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# def nan_sort(mask, x):\n", - "# x = np.sort(np.where(mask, x, -np.inf), axis=1)[:, ::-1]\n", - "# x[x == -np.inf] = np.nan\n", - "# return x\n", - "\n", - "\n", - "\n", - "\n", - "# nan_sort(\n", - "# mask, P\n", - "# )\n", - " \n", - "# # nan_shift2(),\n", - "# # )\n", - "\n", - "# # %timeit nanshift(np.where(mask, P, np.nan))\n", - "# # %timeit nan_shift2()\n", - "\n", - "mpcalc.most_unstable_cape_cin" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/tmp/ipykernel_93894/615816705.py:35: UserWarning: Interpolation point out of data bounds encountered\n", - " ep_, et_, etd_, ptp_ = [x.m for x in mpcalc.parcel_profile_with_lcl(\n" - ] - }, - { - "ename": "AssertionError", - "evalue": "\nNot equal to tolerance rtol=0.001, atol=0\n\nMismatched elements: 1 / 41 (2.44%)\nMax absolute difference: 0.37127686\nMax relative difference: 0.00127804\n x: array([290.5055 , 290.87677, 284.92957, 281.49792, 280.63437, 280.51614,\n 280.40613, 280.20428, 278.98447, 276.37292, 274.80176, 273.24625,\n 271.04562, 268.35522, 265.32175, 262.16504, 258.94226, 256.2096 ,...\n y: array([290.505493, 290.505493, 284.929565, 281.497925, 280.634369,\n 280.516144, 280.406128, 280.204285, 278.984467, 276.372925,\n 274.801758, 273.246246, 271.045624, 268.355225, 265.321747,...", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[9], line 58\u001b[0m\n\u001b[1;32m 55\u001b[0m assert_allclose(et[i], et_, rtol\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1e-3\u001b[39m)\n\u001b[1;32m 57\u001b[0m \u001b[38;5;66;03m# print(etd[i], etd_, sep='\\n')\u001b[39;00m\n\u001b[0;32m---> 58\u001b[0m \u001b[43massert_allclose\u001b[49m\u001b[43m(\u001b[49m\u001b[43metd\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43metd_\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e-3\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 60\u001b[0m \u001b[38;5;66;03m# print(ptp[i], ptp_, sep=\"\\n\")\u001b[39;00m\n\u001b[1;32m 61\u001b[0m assert_allclose(ptp[i], ptp_, rtol\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1e-2\u001b[39m)\n", - " \u001b[0;31m[... skipping hidden 1 frame]\u001b[0m\n", - "File \u001b[0;32m/usr/lib/python3.12/contextlib.py:81\u001b[0m, in \u001b[0;36mContextDecorator.__call__..inner\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m 78\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(func)\n\u001b[1;32m 79\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21minner\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwds):\n\u001b[1;32m 80\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_recreate_cm():\n\u001b[0;32m---> 81\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwds\u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/nzthermo/.venv/lib/python3.12/site-packages/numpy/testing/_private/utils.py:797\u001b[0m, in \u001b[0;36massert_array_compare\u001b[0;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf, strict)\u001b[0m\n\u001b[1;32m 793\u001b[0m err_msg \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;241m.\u001b[39mjoin(remarks)\n\u001b[1;32m 794\u001b[0m msg \u001b[38;5;241m=\u001b[39m build_err_msg([ox, oy], err_msg,\n\u001b[1;32m 795\u001b[0m verbose\u001b[38;5;241m=\u001b[39mverbose, header\u001b[38;5;241m=\u001b[39mheader,\n\u001b[1;32m 796\u001b[0m names\u001b[38;5;241m=\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mx\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124my\u001b[39m\u001b[38;5;124m'\u001b[39m), precision\u001b[38;5;241m=\u001b[39mprecision)\n\u001b[0;32m--> 797\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAssertionError\u001b[39;00m(msg)\n\u001b[1;32m 798\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m:\n\u001b[1;32m 799\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mtraceback\u001b[39;00m\n", - "\u001b[0;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=0.001, atol=0\n\nMismatched elements: 1 / 41 (2.44%)\nMax absolute difference: 0.37127686\nMax relative difference: 0.00127804\n x: array([290.5055 , 290.87677, 284.92957, 281.49792, 280.63437, 280.51614,\n 280.40613, 280.20428, 278.98447, 276.37292, 274.80176, 273.24625,\n 271.04562, 268.35522, 265.32175, 262.16504, 258.94226, 256.2096 ,...\n y: array([290.505493, 290.505493, 284.929565, 281.497925, 280.634369,\n 280.516144, 280.406128, 280.204285, 278.984467, 276.372925,\n 274.801758, 273.246246, 271.045624, 268.355225, 265.321747,..." - ] - } - ], - "source": [ - "(P, T, Td), (Z, Y, X), (latitude, longitude, timestamp, extent) = get_data(\n", - " x = slice(None, None, 45), \n", - " y = slice(None, None, 45),\n", - ")\n", - "\n", - "\n", - "\n", - "\n", - "idx = most_unstable_parcel_index(P, T, Td)\n", - "\n", - "mask = (idx[:, None] <= np.arange(Z)[None, :])\n", - "\n", - "\n", - "\n", - "sort = np.argsort(np.where(mask, P, -np.inf))\n", - "\n", - "\n", - "P_layer, T_layer, Td_layer = np.take_along_axis(\n", - " np.asarray([\n", - " np.where(mask, P, np.nan), \n", - " np.where(mask, T, np.nan), \n", - " np.where(mask, Td, np.nan),\n", - " ]),\n", - " sort[np.newaxis, :, :],\n", - " axis=-1,\n", - ")[:, :, ::-1]\n", - "\n", - "\n", - "# environment pressure, environment temperature, environment dewpoint, parcel temperature profile\n", - "ep, et, etd, ptp = parcel_profile_with_lcl(P_layer, T_layer, Td_layer)\n", - "from numpy.testing import assert_allclose\n", - "# print(ep, et, etd, ptp, sep='\\n')\n", - "for i in range(0, Y*X, 5): \n", - " parcel_idx = idx[i]\n", - " ep_, et_, etd_, ptp_ = [x.m for x in mpcalc.parcel_profile_with_lcl(\n", - " P[parcel_idx:] * Pa,\n", - " T[i, parcel_idx:] * K,\n", - " Td[i, parcel_idx:] * K,\n", - " )]\n", - " # nan_tail = [np.nan]*(Z+1 - len(ep_))\n", - " nan_tail = [np.nan] * (Z+1 - len(ep_))\n", - "\n", - " ep_, et_, etd_, ptp_ = (\n", - " np.array([*x] + nan_tail) for x in (ep_, et_, etd_, ptp_)\n", - " )\n", - "\n", - " etd_[0] = Td[i, parcel_idx]\n", - " ptp_[0] = T[i, parcel_idx]\n", - " ep_[0] = P[parcel_idx]\n", - " et_[0] = T[i, parcel_idx]\n", - "\n", - " # assert_allclose(ep[i], ep_, rtol=1e-3)\n", - "\n", - " # print(et[i], et_, sep='\\n')\n", - " assert_allclose(et[i], et_, rtol=1e-3)\n", - "\n", - " # print(etd[i], etd_, sep='\\n')\n", - " assert_allclose(etd[i], etd_, rtol=1e-3)\n", - "\n", - " # print(ptp[i], ptp_, sep=\"\\n\")\n", - " assert_allclose(ptp[i], ptp_, rtol=1e-2)\n", - "\n", - " # break\n", - " # assert_allclose(etd[i], etd_, rtol=1e-1)\n", - " assert_allclose(ptp[i], ptp_, rtol=1e-2)\n", - " # break\n" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": ".venv", - "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.12.3" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -}