diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index f1a8586d..837c1a06 100644 --- a/docs/internal_api/index.rst +++ b/docs/internal_api/index.rst @@ -53,12 +53,8 @@ Internal Functionality Helpers geocat.comp.meteorology._dewtemp - geocat.comp.meteorology._heat_index - geocat.comp.meteorology._nws_eqn geocat.comp.meteorology._relhum_ice geocat.comp.meteorology._relhum_water - - geocat.comp.meteorology._xheat_index diff --git a/docs/release-notes.rst b/docs/release-notes.rst index 9a250cf4..6e368c39 100644 --- a/docs/release-notes.rst +++ b/docs/release-notes.rst @@ -16,6 +16,10 @@ New Features ^^^^^^^^^^^^ * Update ``calendar_average`` and ``climate_anomaly`` to allow for monthly data with non-uniform spacing by `Katelyn FitzGerald`_ in (:pr:`805`) +Bug Fixes +^^^^^^^^^ +* Rewrite ``heat_index`` so that the calculation is be more in line with the NWS implementation and removes internal functions ``_xheat_index`` and ``_heat_index`` by `Anissa Zacharias`_ in (:pr:`807`) + v2025.12.1 (December 17, 2025) ------------------------------ This release fixes a bug in ``nmse`` and adds a new function, ``delta_pressure_hybrid``. diff --git a/geocat/comp/meteorology.py b/geocat/comp/meteorology.py index 612a8ba2..9a1f4aac 100644 --- a/geocat/comp/meteorology.py +++ b/geocat/comp/meteorology.py @@ -46,123 +46,6 @@ def _dewtemp( return tdk -def _heat_index( - temperature: np.ndarray, - relative_humidity: typing.Union[np.ndarray, xr.DataArray, list, float], - alternate_coeffs: bool = False, -) -> np.ndarray: - """Compute the 'heat index' as calculated by the National Weather Service. - - Internal function for heat_index - - Parameters - ---------- - temperature : ndarray - temperature(s) in Fahrenheit - - relative_humidity : ndarray, :class:`xarray.DataArray`, :class:`list`, :class:`float` - relative humidity as a percentage. Must be the same shape as - ``temperature`` - - alternate_coeffs : bool, optional - flag to use alternate set of coefficients appropriate for - temperatures from 70F to 115F and humidities between 0% and 80% - - Returns - ------- - heatindex : ndarray - Calculated heat index. Same shape as ``temperature`` - - - See Also - -------- - Related GeoCAT Functions: - :func:`heat_index`, - :func:`_xheat_index` - - Related NCL Functions: - `heat_index_nws `__ - """ - # Default coefficients for (t>=80F) and (40 crit[0], - _nws_eqn(coeffs, temperature, relative_humidity), - heatindex, - ) - - # adjustments - heatindex = xr.where( - np.logical_and( - relative_humidity < 13, - np.logical_and(temperature > 80, temperature < 112), - ), - heatindex - - ((13 - relative_humidity) / 4) - * np.sqrt((17 - abs(temperature - 95)) / 17), - heatindex, - ) - - heatindex = xr.where( - np.logical_and( - relative_humidity > 85, - np.logical_and(temperature > 80, temperature < 87), - ), - heatindex - + ((relative_humidity - 85.0) / 10.0) * ((87.0 - temperature) / 5.0), - heatindex, - ) - - return heatindex - - def _nws_eqn(coeffs, temp, rel_hum): """Helper function to compute the heat index. @@ -189,8 +72,7 @@ def _nws_eqn(coeffs, temp, rel_hum): See Also -------- Related GeoCAT Functions: - :func:`heat_index`, - :func:`_heat_index` + :func:`heat_index` Related NCL Functions: `heat_index_nws `__, @@ -367,128 +249,6 @@ def _relhum_water( return rh -def _xheat_index( - temperature: xr.DataArray, - relative_humidity: xr.DataArray, - alternate_coeffs: bool = False, -) -> tuple([xr.DataArray, int]): - """Compute the 'heat index' as calculated by the National Weather Service. - - Internal function for heat_index for dask - - Parameters - ---------- - temperature : :class:`xarray.DataArray` - temperature(s) in Fahrenheit - - relative_humidity : :class:`xarray.DataArray` - relative humidity as a percentage. Must be the same shape as - ``temperature`` - - alternate_coeffs : bool, optional - flag to use alternate set of coefficients appropriate for - temperatures from 70F to 115F and humidities between 0% and 80% - - Returns - ------- - heatindex : :class:`xarray.DataArray` - Calculated heat index. Same shape as ``temperature`` - - eqtype : :class:`int` - version of equations used, for xarray attrs output - - See Also - -------- - Related GeoCAT Functions: - :func:`heat_index` - :func:`_heat_index` - - Related NCL Functions: - `heat_index_nws `__, - """ - # Default coefficients for (t>=80F) and (40 crit[0], - _nws_eqn(coeffs, temperature, relative_humidity), - heatindex, - ) - - # adjustments - heatindex = xr.where( - np.logical_and( - relative_humidity < 13, - np.logical_and(temperature > 80, temperature < 112), - ), - heatindex - - ((13 - relative_humidity) / 4) - * np.sqrt((17 - abs(temperature - 95)) / 17), - heatindex, - ) - - heatindex = xr.where( - np.logical_and( - relative_humidity > 85, - np.logical_and(temperature > 80, temperature < 87), - ), - heatindex - + ((relative_humidity - 85.0) / 10.0) * ((87.0 - temperature) / 5.0), - heatindex, - ) - - return heatindex, eqtype - - def dewtemp( temperature: typing.Union[np.ndarray, xr.DataArray, list, float], relative_humidity: typing.Union[np.ndarray, xr.DataArray, list, float], @@ -569,7 +329,6 @@ def heat_index( The computation of the heat index is a refinement of a result obtained by multiple regression analysis carried out by Lans P. Rothfusz and described in a 1990 National Weather Service (NWS) Technical Attachment (SR 90-23). - All values less that 40F/4.4C/277.65K are set to the ambient temperature. In practice, the Steadman formula is computed first and the result averaged with the temperature. If this heat index value is 80 degrees F or higher, the full regression equation along with any adjustment as described above @@ -607,74 +366,126 @@ def heat_index( See Also -------- - Related GeoCAT Functions: - :func:`_heat_index` - :func:`_xheat_index` - Related NCL Functions: `heat_index_nws `__ """ - inputs = [temperature, relative_humidity] - # ensure all inputs same size - if not (np.shape(x) == np.shape(inputs[0]) for x in inputs): - raise ValueError("heat_index: dimensions of inputs are not the same") - - # Get input types + # get input types in_types = [type(item) for item in inputs] - # run dask compatible version if input is xarray - if xr.DataArray in in_types: - # check all inputs are xarray.DataArray - if not all(x == xr.DataArray for x in in_types): - raise TypeError("heat_index: if using xarray, all inputs must be xarray") + # ensure all same input type + if len(set(in_types)) != 1: + raise TypeError( + f"heat_index: input types are not the same, received {in_types}" + ) - # input validation on relative humidity - if any(relative_humidity.data.ravel() < 0) or any( - relative_humidity.data.ravel() > 100 - ): - raise ValueError('heat_index: invalid values for relative humidity') + # if inputs not xarray or numpy (float, int, list), elevate to np + if not isinstance(temperature, xr.DataArray) and not isinstance( + temperature, np.ndarray + ): + try: + temperature = np.asarray(temperature) + relative_humidity = np.asarray(relative_humidity) + except [ValueError, TypeError] as e: + raise TypeError(f"heat_index: cannot convert input to numpy array, {e}") - # Check if relative humidity fractional - if all(relative_humidity.data.ravel() < 1): - warnings.warn("WARNING: rh must be %, not fractional; All rh are < 1") + # ensure all inputs same size + if not temperature.shape == relative_humidity.shape: + raise ValueError( + f"heat_index: dimensions of inputs are not the same, received {temperature.shape}, {relative_humidity.shape}" + ) - # call internal computation function - heatindex, eqtype = _xheat_index( - temperature, relative_humidity, alternate_coeffs + # check relative humidity within valid range 0 to 100 + if relative_humidity.min() < 0 or relative_humidity.max() > 100: + raise ValueError( + 'heat_index: invalid values for relative humidity, expected all between 0 and 100' ) - # set xarray attributes + # check if relative humidity fractional + if relative_humidity.max() <= 1: + warnings.warn("WARNING: rh must be %, not fractional; All rh are <= 1") + + # Default coefficients for (t>=80F) and (40 79.0, + # but the equation page says "80 degrees F or higher", we're going to + # use Rothfusz for > 79 F for default coeffs to match the published table + # https://www.weather.gov/safety/heat-tools + # https://www.weather.gov/images/safety/heatindexchart-650.jpg + heatindex = xr.where(heatindex >= crit, rothfusz, heatindex) + + # adjustment for rh <= 13, 80F <= t <= 112F + heatindex = xr.where( + np.logical_and( + relative_humidity <= 13, + np.logical_and(temperature >= 80, temperature <= 112), + ), + heatindex + - ((13 - relative_humidity) / 4) + * np.sqrt(abs((17 - abs(temperature - 95))) / 17), + heatindex, + ) + # adjustment for rh > 85 and 80F <= t <= 87F + heatindex = xr.where( + np.logical_and( + relative_humidity > 85, + np.logical_and(temperature >= 80, temperature <= 87), + ), + heatindex + ((relative_humidity - 85.0) / 10.0) * ((87.0 - temperature) / 5.0), + heatindex, + ) + + # if elevated to np array, return to original type + if in_types[0] in {list, float, int} and isinstance(heatindex, np.ndarray): + heatindex = heatindex.tolist() + + # if xarray, set attributes + if in_types[0] == xr.DataArray: heatindex.attrs['long_name'] = "heat index: NWS" heatindex.attrs['units'] = "F" heatindex.attrs['www'] = ( "https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml" ) heatindex.attrs['info'] = "appropriate for shady locations with no wind" - - if eqtype == 1: - heatindex.attrs['tag'] = ( - "NCL: heat_index_nws; (Steadman+t)*0.5 and Rothfusz" - ) - else: - heatindex.attrs['tag'] = "NCL: heat_index_nws; (Steadman+t)*0.5" - - else: - # ensure in numpy array for function call - temperature = np.atleast_1d(temperature) - relative_humidity = np.atleast_1d(relative_humidity) - - # input validation on relative humidity - if any(relative_humidity.ravel() < 0) or any(relative_humidity.ravel() > 100): - raise ValueError('heat_index: invalid values for relative humidity') - - # Check if relative humidity fractional - if all(relative_humidity.ravel() < 1): - warnings.warn("WARNING: rh must be %, not fractional; All rh are < 1") - - # function call for non-dask/xarray - heatindex = _heat_index(temperature, relative_humidity, alternate_coeffs) + heatindex.attrs['tag'] = "heat_index: (Steadman+t)*0.5 and Rothfusz" return heatindex diff --git a/test/heat-index.csv b/test/heat-index.csv new file mode 100644 index 00000000..e420068d --- /dev/null +++ b/test/heat-index.csv @@ -0,0 +1,16 @@ +79.9,80.3,80.8,81.3,81.8,82.4,83,83.6,84.2,84.9,86.3,87.8,89.3 +81.5,82.1,82.8,83.6,84.4,85.4,86.4,87.6,88.8,90.1,92,94,96 +83.3,84.1,85.1,86.3,87.5,88.9,90.5,92.2,94,95.9,98.3,100.9,103.6 +85.4,86.6,87.9,89.4,91.1,93,95.1,97.3,99.8,102.5,105.4,108.5,111.8 +87.9,89.4,91,93,95.1,97.6,100.2,103.1,106.3,109.6,113.2,117.1,121.2 +90.7,92.5,94.6,97,99.7,102.7,105.9,109.5,113.3,117.5,121.9,126.6,131.6 +93.8,96,98.5,101.4,104.7,108.3,112.2,116.4,121,126,131.3,136.9,142.8 +97.2,99.8,102.9,106.3,110.2,114.4,119,124,129.4,135.2,141.3,147.9,154.8 +100.9,104,107.6,111.7,116.1,121,126.4,132.1,138.3,145,152.1,159.6,167.5 +104.9,108.6,112.8,117.4,122.6,128.2,134.3,140.9,147.9,155.5,163.5,172,181 +109.3,113.5,118.3,123.6,129.5,135.9,142.8,150.2,158.2,166.7,175.7,185.2,195.3 +113.9,118.8,124.3,130.3,136.9,144.1,151.8,160.1,169,178.5,188.5,199.2,210.4 +118.9,124.4,130.6,137.4,144.8,152.8,161.4,170.7,180.5,191,202.1,213.8,226.2 +124.2,130.4,137.3,144.9,153.1,162,171.6,181.8,192.6,204.2,216.4,229.2,242.7 +129.8,136.7,144.4,152.8,161.9,171.7,182.3,193.5,205.4,218,231.3,245.4,260.1 +135.7,143.4,152,161.2,171.2,182,193.5,205.8,218.8,232.5,247,262.2,278.2 diff --git a/test/test_dask.py b/test/test_dask.py index 4fd19333..98c0a91b 100644 --- a/test/test_dask.py +++ b/test/test_dask.py @@ -102,27 +102,26 @@ def test_dewtemp_dask(self): assert np.allclose(out - 273.15, dt_2, atol=0.1) def test_heat_index_dask(self): - ncl_gt_1 = [ - 137.36142, - 135.86795, - 104.684456, - 131.25621, - 105.39449, - 79.78999, - 83.57511, - 59.965, - 30.0, + hi_ncl_alt = [ + 76.13114, + 75.12854, + 99.43573, + 104.93261, + 93.73293, + 104.328705, + 123.23398, + 150.34001, + 106.87023, ] - - t1 = np.array([104, 100, 92, 92, 86, 80, 80, 60, 30]) - rh1 = np.array([55, 65, 60, 90, 90, 40, 75, 90, 50]) + t1 = np.array([75, 80, 85, 90, 95, 100, 105, 110, 115]) + rh1 = np.array([75, 15, 80, 65, 25, 30, 40, 50, 5]) t = xr.DataArray(t1).chunk(3) rh = xr.DataArray(rh1).chunk(3) - out = heat_index(t, rh) + out = heat_index(t, rh, alternate_coeffs=True) assert isinstance(out.data, dask.array.Array) - assert np.allclose(out, ncl_gt_1, atol=0.005) + assert np.allclose(out, hi_ncl_alt, atol=0.005) def test_relhum_dask(self): p_def = [ diff --git a/test/test_meteorology.py b/test/test_meteorology.py index 1bf4fa12..74a7a030 100644 --- a/test/test_meteorology.py +++ b/test/test_meteorology.py @@ -119,19 +119,8 @@ def test_xarray_type_error(self) -> None: class Test_heat_index: # set up ground truths - ncl_gt_1 = [ - 137.36142, - 135.86795, - 104.684456, - 131.25621, - 105.39449, - 79.78999, - 83.57511, - 59.965, - 30.0, - ] - ncl_gt_2 = [ - 68.585, + # use ncl for alt coefficient testing within unchanged ranges + hi_ncl_alt = [ 76.13114, 75.12854, 99.43573, @@ -142,63 +131,60 @@ class Test_heat_index: 150.34001, 106.87023, ] + t_alt = np.array([75, 80, 85, 90, 95, 100, 105, 110, 115]) + rh_alt = np.array([75, 15, 80, 65, 25, 30, 40, 50, 5]) - t1 = np.array([104, 100, 92, 92, 86, 80, 80, 60, 30]) - rh1 = np.array([55, 65, 60, 90, 90, 40, 75, 90, 50]) + t_nws = [80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110] + rh_nws = [40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100] + t_nws, rh_nws = np.meshgrid(t_nws, rh_nws) - t2 = np.array([70, 75, 80, 85, 90, 95, 100, 105, 110, 115]) - rh2 = np.array([10, 75, 15, 80, 65, 25, 30, 40, 50, 5]) + # used javascript to query NWS form + @pytest.fixture(scope="class") + def hi_nws(self): + try: + return np.genfromtxt("heat-index.csv", delimiter=",").T + except Exception: + return np.genfromtxt("test/heat-index.csv", delimiter=",").T - def test_numpy_input(self) -> None: - assert np.allclose( - heat_index(self.t1, self.rh1, False), self.ncl_gt_1, atol=0.005 - ) + def test_numpy_input(self, hi_nws) -> None: + hi = heat_index(self.t_nws, self.rh_nws, False) + assert np.allclose(hi.round(1), hi_nws) def test_multi_dimensional_input(self) -> None: assert np.allclose( - heat_index(self.t2.reshape(2, 5), self.rh2.reshape(2, 5), True), - np.asarray(self.ncl_gt_2).reshape(2, 5), + heat_index(self.t_alt.reshape(3, 3), self.rh_alt.reshape(3, 3), True), + np.asarray(self.hi_ncl_alt).reshape(3, 3), atol=0.005, ) def test_alt_coef(self) -> None: assert np.allclose( - heat_index(self.t2, self.rh2, True), self.ncl_gt_2, atol=0.005 + heat_index(self.t_alt, self.rh_alt, True), self.hi_ncl_alt, atol=0.005 ) def test_xarray_alt_coef(self) -> None: - assert np.allclose( - heat_index(xr.DataArray(self.t2), xr.DataArray(self.rh2), True), - self.ncl_gt_2, - atol=0.005, - ) + hi_alt = heat_index(xr.DataArray(self.t_alt), xr.DataArray(self.rh_alt), True) + assert np.allclose(hi_alt, self.hi_ncl_alt, atol=0.005) def test_float_input(self) -> None: assert np.allclose(heat_index(80, 75), 83.5751, atol=0.005) - def test_list_input(self) -> None: - assert np.allclose( - heat_index(self.t1.tolist(), self.rh1.tolist()), self.ncl_gt_1, atol=0.005 - ) + def test_list_input(self, hi_nws) -> None: + hi = heat_index(self.t_nws.tolist(), self.rh_nws.tolist()) + hi = [[round(x, 1) for x in r] for r in hi] + assert np.allclose(hi, hi_nws) - def test_xarray_input(self) -> None: - t = xr.DataArray(self.t1) - rh = xr.DataArray(self.rh1) + def test_xarray_input(self, hi_nws) -> None: + t = xr.DataArray(self.t_nws) + rh = xr.DataArray(self.rh_nws) - assert np.allclose(heat_index(t, rh), self.ncl_gt_1, atol=0.005) - - def test_alternate_xarray_tag(self) -> None: - t = xr.DataArray([15, 20]) - rh = xr.DataArray([15, 20]) - - out = heat_index(t, rh) - assert out.tag == "NCL: heat_index_nws; (Steadman+t)*0.5" + assert np.allclose(heat_index(t, rh).round(1), hi_nws) def test_rh_warning(self) -> None: with pytest.warns(UserWarning): heat_index([50, 80, 90], [0.1, 0.2, 0.5]) - def test_rh_valid(self) -> None: + def test_rh_invalid(self) -> None: with pytest.raises(ValueError): heat_index([50, 80, 90], [-1, 101, 50]) @@ -212,11 +198,15 @@ def test_xarray_rh_valid(self) -> None: def test_xarray_type_error(self) -> None: with pytest.raises(TypeError): - heat_index(self.t1, xr.DataArray(self.rh1)) + heat_index(self.t_nws, xr.DataArray(self.rh_nws)) def test_dims_error(self) -> None: with pytest.raises(ValueError): - heat_index(self.t1[:10], self.rh1[:8]) + heat_index(self.t_nws[:10], self.rh_nws[:8]) + + def test_bad_list_input(self): + with pytest.raises(ValueError): + heat_index(np.asarray([[85, 85], [85]]), np.asarray([[60, 60], [60]])) class Test_relhum: