diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst index 2bc99905..036bbbf5 100644 --- a/docs/internal_api/index.rst +++ b/docs/internal_api/index.rst @@ -43,6 +43,8 @@ Internal API geocat.comp.interpolation._post_interp_multidim + geocat.comp.meteorology._delta_pressure1D + geocat.comp.meteorology._dewtemp geocat.comp.meteorology._heat_index diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst index a6da25f2..5edbb773 100644 --- a/docs/user_api/index.rst +++ b/docs/user_api/index.rst @@ -58,6 +58,7 @@ Meteorology :toctree: ./generated/ actual_saturation_vapor_pressure + delta_pressure dewtemp heat_index max_daylight @@ -98,7 +99,6 @@ GeoCAT-comp routines from GeoCAT-f2py :nosignatures: :toctree: ./generated/ - dpres_plevel grid_to_triple linint1 linint2 @@ -109,6 +109,16 @@ GeoCAT-comp routines from GeoCAT-f2py rgrid2rcm triple_to_grid +NCL Function Name Wrappers +-------------------------- +.. currentmodule:: geocat.comp +.. autosummary:: + :nosignatures: + :toctree: ./generated/ + + meteorology.dpres_plev + + Deprecated Functions -------------------- Climatologies diff --git a/src/geocat/comp/__init__.py b/src/geocat/comp/__init__.py index 8c29dd72..139b0f78 100644 --- a/src/geocat/comp/__init__.py +++ b/src/geocat/comp/__init__.py @@ -8,7 +8,7 @@ from .meteorology import (dewtemp, heat_index, relhum, relhum_ice, relhum_water, actual_saturation_vapor_pressure, max_daylight, psychrometric_constant, saturation_vapor_pressure, - saturation_vapor_pressure_slope) + saturation_vapor_pressure_slope, delta_pressure) from .spherical import decomposition, recomposition, scale_voronoi from .stats import eofunc, eofunc_eofs, eofunc_pcs, eofunc_ts, pearson_r # bring all functions from geocat.f2py into the geocat.comp namespace diff --git a/src/geocat/comp/meteorology.py b/src/geocat/comp/meteorology.py index e8065fe7..1b73af08 100644 --- a/src/geocat/comp/meteorology.py +++ b/src/geocat/comp/meteorology.py @@ -1379,3 +1379,142 @@ def saturation_vapor_pressure_slope( tfill) return svp_slope + + +def _delta_pressure1D(pressure_lev, surface_pressure): + """Helper function for `delta_pressure`. Calculates the pressure layer + thickness (delta pressure) of a one-dimensional pressure level array. + + Returns an array of length matching `pressure_lev`. + + Parameters + ---------- + pressure_lev : :class:`numpy.ndarray` + The pressure level array. May be in ascending or descending order. + Must have the same units as `surface_pressure`. + + surface_pressure : :class:`float` + The scalar surface pressure. Must have the same units as + `pressure_lev`. + + Returns + ------- + delta_pressure : :class:`numpy.ndarray` + The pressure layer thickness array. Shares dimensions and units of + `pressure_lev`. + """ + pressure_top = min(pressure_lev) + + # Safety checks + if pressure_top <= 0: + warnings.warn("'pressure_lev` values must all be positive.") + if pressure_top > surface_pressure: + warnings.warn( + "`surface_pressure` must be greater than minimum `pressure_lev` value." + ) + + # Sort so pressure increases (array goes from top of atmosphere to bottom) + is_pressuredecreasing = pressure_lev[1] < pressure_lev[0] + if is_pressuredecreasing: + pressure_lev = np.flip(pressure_lev) + + # Calculate delta pressure + delta_pressure = np.empty_like(pressure_lev) + + delta_pressure[0] = (pressure_lev[0] + + pressure_lev[1]) / 2 - pressure_top # top level + delta_pressure[1:-1] = [ + (a - b) / 2 for a, b in zip(pressure_lev[2:], pressure_lev[:-1]) + ] + delta_pressure[-1] = surface_pressure - ( + pressure_lev[-1] + pressure_lev[-2]) / 2 # bottom level + + # Return delta_pressure to original order + if is_pressuredecreasing: + delta_pressure = np.flip(delta_pressure) + + return delta_pressure + + +def delta_pressure(pressure_lev, surface_pressure): + """Calculates the pressure layer thickness (delta pressure) of a constant + pressure level coordinate system. + + Returns an array of shape matching (`surface_pressure`, `pressure_lev`). + + Parameters + ---------- + pressure_lev : :class:`numpy.ndarray`, :class:'xarray.DataArray` + The pressure level array. May be in ascending or descending order. + Must have the same units as `surface_pressure`. + surface_pressure : :class:`np.Array`, :class:'xr.DataArray` + The scalar or N-dimensional surface pressure array. Must have the same + units as `pressure_lev`. + + Returns + ------- + delta_pressure : :class:`numpy.ndarray`, :class:'xarray.DataArray` + The pressure layer thickness array. Shares units with `pressure_lev`. + If `surface_pressure` is scalar, shares dimensions with + `pressure_level`. If `surface_pressure` is an array than the returned + array will have an additional dimension [e.g. (lat, lon, time) becomes + (lat, lon, time, lev)]. + + See Also + -------- + Related NCL Functions: + `dpres_plev `__ + """ + # Get original array types + type_surface_pressure = type( + surface_pressure + ) # save type for delta_pressure to same type as surface_pressure at end + type_pressure_level = type(pressure_lev) + + # Preserve attributes for Xarray + if type_surface_pressure == xr.DataArray: + da_coords = dict(surface_pressure.coords) + da_attrs = dict(surface_pressure.attrs) + da_dims = surface_pressure.dims + if type_pressure_level == xr.DataArray: + da_attrs = dict( + pressure_lev.attrs) # Overwrite attributes to match pressure_lev + + # Calculate delta pressure + if np.isscalar(surface_pressure): # scalar case + delta_pressure = _delta_pressure1D(pressure_lev, surface_pressure) + else: # multi-dimensional cases + shape = surface_pressure.shape + delta_pressure_shape = shape + (len(pressure_lev), + ) # preserve shape for reshaping + + surface_pressure_flattened = np.ravel( + surface_pressure) # flatten to avoid nested for loops + delta_pressure = [ + _delta_pressure1D(pressure_lev, e) + for e in surface_pressure_flattened + ] + + delta_pressure = np.array(delta_pressure).reshape(delta_pressure_shape) + + # If passed in an Xarray array, return an Xarray array + # Change this to return a dataset that has both surface pressure and delta pressure? + if type_surface_pressure == xr.DataArray: + da_coords['lev'] = pressure_lev.values + da_dims = da_dims + ("lev",) + da_attrs.update({"long name": "pressure layer thickness"}) + delta_pressure = xr.DataArray(delta_pressure, + coords=da_coords, + dims=da_dims, + attrs=da_attrs, + name="delta pressure") + + return delta_pressure + + +def dpres_plev(pressure_lev, surface_pressure): + return delta_pressure(pressure_lev, surface_pressure) + + +_dpres_plev_doc_str = f".. attention:: This method is a wrapper for `delta_pressure `_.\n\n {delta_pressure.__doc__}" +setattr(dpres_plev, '__doc__', _dpres_plev_doc_str) diff --git a/test/test_meteorology.py b/test/test_meteorology.py index 6149ec9a..22f5a484 100644 --- a/test/test_meteorology.py +++ b/test/test_meteorology.py @@ -16,14 +16,14 @@ from src.geocat.comp.meteorology import ( dewtemp, heat_index, relhum, relhum_ice, relhum_water, actual_saturation_vapor_pressure, max_daylight, psychrometric_constant, - saturation_vapor_pressure, saturation_vapor_pressure_slope) + saturation_vapor_pressure, saturation_vapor_pressure_slope, + delta_pressure) else: - from geocat.comp.meteorology import (dewtemp, heat_index, relhum, - relhum_ice, relhum_water, - actual_saturation_vapor_pressure, - max_daylight, psychrometric_constant, - saturation_vapor_pressure, - saturation_vapor_pressure_slope) + from geocat.comp.meteorology import ( + dewtemp, heat_index, relhum, relhum_ice, relhum_water, + actual_saturation_vapor_pressure, max_daylight, psychrometric_constant, + saturation_vapor_pressure, saturation_vapor_pressure_slope, + delta_pressure) class Test_dewtemp(unittest.TestCase): @@ -645,3 +645,89 @@ def test_dask_lazy(self): assert isinstance((saturation_vapor_pressure_slope(tempf)).data, dask.array.Array) + + +class TestDeltaPressure(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.pressure_lev = np.array([1, 5, 100, 1000]) + cls.pressure_lev_da = xr.DataArray(cls.pressure_lev) + cls.pressure_lev_da.attrs = { + "long name": "pressure level", + "units": "hPa", + "direction": "descending" + } + + cls.surface_pressure_scalar = 1018 + cls.surface_pressure_1D = np.array([1018, 1019]) + cls.surface_pressure_2D = np.array([[1018, 1019], [1017, 1019.5]]) + cls.surface_pressure_3D = np.array([[[1018, 1019], [1017, 1019.5]], + [[1019, 1020], [1018, 1020.5]]]) + + coords = {'time': [1, 2], 'lat': [3, 4], 'lon': [5, 6]} + dims = ["time", "lat", "lon"] + attrs = {"long name": "surface pressure", "units": "hPa"} + cls.surface_pressure_3D_da = xr.DataArray(cls.surface_pressure_3D, + coords=coords, + dims=dims, + attrs=attrs) + + def test_delta_pressure1D(self): + pressure_lev = [float(i) for i in self.pressure_lev] + pressure_top = min(pressure_lev) + delta_p = delta_pressure(pressure_lev, self.surface_pressure_scalar) + self.assertEqual(sum(delta_p), + self.surface_pressure_scalar - pressure_top) + + def test_negative_pressure_warning(self): + pressure_lev_negative = self.pressure_lev.copy() + pressure_lev_negative[0] = -5 + with self.assertWarns(Warning): + delta_p = delta_pressure(pressure_lev_negative, + self.surface_pressure_scalar) + + def test_relative_pressure_warning(self): + surface_pressure_low = 0.5 + with self.assertWarns(Warning): + delta_p = delta_pressure(self.pressure_lev, surface_pressure_low) + + def test_output_type(self): + delta_pressure_da = delta_pressure(self.pressure_lev_da, + self.surface_pressure_3D_da) + self.assertIsInstance(delta_pressure_da, xr.DataArray) + + delta_pressure_np = delta_pressure(self.pressure_lev, + self.surface_pressure_3D) + self.assertIsInstance(delta_pressure_np, np.ndarray) + + def test_output_dimensions(self): + delta_pressure_scalar = delta_pressure(self.pressure_lev, + self.surface_pressure_scalar) + self.assertEqual(delta_pressure_scalar.shape, (4,)) + + delta_pressure_1D = delta_pressure(self.pressure_lev, + self.surface_pressure_1D) + self.assertEqual(delta_pressure_1D.shape, (2, 4)) + + delta_pressure_2D = delta_pressure(self.pressure_lev, + self.surface_pressure_2D) + self.assertEqual(delta_pressure_2D.shape, (2, 2, 4)) + + delta_pressure_3D = delta_pressure(self.pressure_lev, + self.surface_pressure_3D) + self.assertEqual(delta_pressure_3D.shape, (2, 2, 2, 4)) + + def test_output_attrs(self): + delta_pressure_da = delta_pressure(self.pressure_lev_da, + self.surface_pressure_3D_da) + for item in self.pressure_lev_da.attrs: + self.assertIn(item, delta_pressure_da.attrs) + + def test_output_coords(self): + delta_pressure_da = delta_pressure(self.pressure_lev_da, + self.surface_pressure_3D_da) + for item in self.surface_pressure_3D_da.coords: + self.assertIn(item, delta_pressure_da.coords) + for item in self.pressure_lev_da.coords: + self.assertIn(item, delta_pressure_da.coords)