From a99c908952722bf2bfc38c75b7ee1b6178e4a505 Mon Sep 17 00:00:00 2001 From: Thomas Rieutord Date: Fri, 19 Jul 2024 12:14:20 +0100 Subject: [PATCH 1/5] Add a boundary layer module with classical ways to estimate boundary layer height --- src/metpy/calc/boundarylayer.py | 305 ++++++++++++++++++++++++++++++++ tests/test_boundarylayer.py | 145 +++++++++++++++ 2 files changed, 450 insertions(+) create mode 100644 src/metpy/calc/boundarylayer.py create mode 100644 tests/test_boundarylayer.py diff --git a/src/metpy/calc/boundarylayer.py b/src/metpy/calc/boundarylayer.py new file mode 100644 index 0000000000..1d43ea0ebe --- /dev/null +++ b/src/metpy/calc/boundarylayer.py @@ -0,0 +1,305 @@ +#!/usr/bin/python +# -*-coding:utf-8 -*- +""" +Contains a collection of boundary layer height estimations. + + +References: +----------- + +[Col14]: Collaud Coen, M., Praz, C., Haefele, A., Ruffieux, D., Kaufmann, P., and Calpini, B. (2014): + Determination and climatology of the planetary boundary layer height above the Swiss plateau by in situ and remote sensing measurements as well as by the COSMO-2 model + Atmos. Chem. Phys., 14, 13205–13221. + +[HL06]: Hennemuth, B., & Lammert, A. (2006): + Determination of the atmospheric boundary layer height from radiosonde and lidar backscatter. + Boundary-Layer Meteorology, 120(1), 181-200. + +[Guo16]: Guo, J., Miao, Y., Zhang, Y., Liu, H., Li, Z., Zhang, W., ... & Zhai, P. (2016): + The climatology of planetary boundary layer height in China derived from radiosonde and reanalysis data. + Atmos. Chem. Phys, 16(20), 13309-13319. + +[Sei00]: Seidel, D. J., Ao, C. O., & Li, K. (2010): + Estimating climatological planetary boundary layer heights from radiosonde observations: Comparison of methods and uncertainty analysis. + Journal of Geophysical Research: Atmospheres, 115(D16). + +[VH96]: Vogelezang, D. H. P., & Holtslag, A. A. M. (1996): + Evaluation and model impacts of alternative boundary-layer height formulations. + Boundary-Layer Meteorology, 81(3-4), 245-269. +""" + +import numpy as np +from copy import deepcopy + +import metpy.calc as mpcalc +import metpy.constants as mpconsts +from metpy.units import units + + +def smooth(val, span): + """Function that calculates the moving average with a given span. + The span is given in number of points on which the average is made. + + Parameters + ---------- + val: array-like + Array of values + span: int + Span of the moving average. The higher the smoother + + Returns + ------- + smoothed_val: array-like + Array of smoothed values + """ + N = len(val) + smoothed_val = deepcopy(val) + for i in range(N): + smoothed_val[i] = np.nanmean(val[i - min(span, i) : i + min(span, N - i)]) + + return smoothed_val + + +def bulk_richardson_number( + height, + potential_temperature, + u, + v, + idxfoot: int = 0, + ustar=0 * units.meter_per_second, +): + r"""Calculate the bulk Richardson number. + + See [VH96], eq. (3): + + .. math:: Ri = (g/\theta) * \frac{(\Delta z)(\Delta \theta)} + {\left(\Delta u)^2 + (\Delta v)^2 + b(u_*)^2} + + Parameters + ---------- + height : `pint.Quantity` + Altitude (metres above ground) of the points in the profile + potential_temperature : `pint.Quantity` + Potential temperature profile + u : `pint.Quantity` + Zonal wind profile + v : `pint.Quantity` + Meridional wind profile + idxfoot : int, optional + The index of the foot point (first trusted measure), defaults to 0. + + Returns + ------- + `pint.Quantity` + Bulk Richardson number profile + """ + + u[0] = 0 * units.meter_per_second + v[0] = 0 * units.meter_per_second + + Dtheta = potential_temperature - potential_temperature[idxfoot] + Du = u - u[idxfoot] + Dv = v - v[idxfoot] + Dz = height - height[idxfoot] + + idx0 = Du**2 + Dv**2 + ustar**2 == 0 + if idx0.sum() > 0: + bRi = np.ones_like(Dtheta) * np.nan * units.dimensionless + bRi[~idx0] = ( + (mpconsts.g / potential_temperature[~idx0]) + * (Dtheta[~idx0] * Dz[~idx0]) + / (Du[~idx0] ** 2 + Dv[~idx0] ** 2 + ustar**2) + ) + else: + bRi = ( + (mpconsts.g / potential_temperature) + * (Dtheta * Dz) + / (Du**2 + Dv**2 + ustar**2) + ) + + return bRi + + +def blh_from_richardson_bulk( + height, + potential_temperature, + u, + v, + smoothingspan: int = 10, + idxfoot: int = 0, + bri_threshold=0.25 * units.dimensionless, + ustar=0.1 * units.meter_per_second, +): + """Calculate atmospheric boundary layer height with the method of + bulk Richardson number. + + It is the height where the bulk Richardson number exceeds a given threshold. + See [VH96, Sei00, Col14, Guo16]. + + Parameters + ---------- + height : `pint.Quantity` + Altitude (metres above ground) of the points in the profile + potential_temperature : `pint.Quantity` + Potential temperature profile + u : `pint.Quantity` + Zonal wind profile + v : `pint.Quantity` + Meridional wind profile + smoothingspan : int, optional + The amount of smoothing (number of points in moving average) + idxfoot : int, optional + The index of the foot point (first trusted measure), defaults to 0. + bri_threshold : `pint.Quantity`, optional + Threshold to exceed to get boundary layer top. Defaults to 0.25 + ustar : `pint.Quantity`, optional + Additional friction term in [VH96]. Defaluts to 0. + + Returns + ------- + blh : `pint.Quantity` + Boundary layer height estimation + """ + bRi = bulk_richardson_number( + height, + smooth(potential_temperature, smoothingspan), + smooth(u, smoothingspan), + smooth(v, smoothingspan), + idxfoot=idxfoot, + ustar=ustar, + ) + + height = height[~np.isnan(bRi)] + bRi = bRi[~np.isnan(bRi)] + + if any(bRi > bri_threshold): + iblh = np.where(bRi > bri_threshold)[0][0] + blh = height[iblh] + else: + blh = np.nan * units.meter + + return blh + + +def blh_from_parcel( + height, + potential_temperature, + smoothingspan: int = 5, + theta0=None, +): + """Calculate atmospheric boundary layer height with the parcel method. + + It is the height where the potential temperature profile reaches its + foot value. + See [Sei00, HL06, Col14]. + + Parameters + ---------- + height : `pint.Quantity` + Altitude (metres above ground) of the points in the profile + potential_temperature : `pint.Quantity` + Potential temperature profile + smoothingspan : int, optional + The amount of smoothing (number of points in moving average) + theta0 : `pint.Quantity`, optional + Value of theta at the foot point (skip unstruted points or add extra term). If not provided, theta[0] is taken. + + Returns + ------- + blh : `pint.Quantity` + Boundary layer height estimation + """ + potential_temperature = smooth(potential_temperature, smoothingspan) + + if theta0 is None: + theta0 = potential_temperature[0] + + if any(potential_temperature > theta0): + iblh = np.where(potential_temperature > theta0)[0][0] + blh = height[iblh] + else: + blh = np.nan * units.meter + + return blh + + +def blh_from_humidity_gradient( + height, + humidity, + smoothingspan: int = 5, + idxfoot: int = 0, +): + """Calculate atmospheric boundary layer height from the relative + humidity gradient + + It is the height where the relative humidity or specific humidity gradient reaches a minimum. + See [Sei00, HL06, Col14]. + + Parameters + ---------- + height : `pint.Quantity` + Altitude (metres above ground) of the points in the profile + humidity : `pint.Quantity` + Humidity (relative or specific) profile + smoothingspan : int, optional + The amount of smoothing (number of points in moving average) + idxfoot : int, optional + The index of the foot point (first trusted measure), defaults to 0. + + Returns + ------- + blh : `pint.Quantity` + Boundary layer height estimation + """ + + dRHdz = mpcalc.first_derivative(smooth(humidity, smoothingspan), x=height) + + dRHdz = dRHdz[idxfoot:] + height = height[idxfoot:] + + iblh = np.argmin(dRHdz) + + return height[iblh] + + +def blh_from_temperature_inversion( + height, + temperature, + smoothingspan: int = 5, + idxfoot: int = 0, +): + """Calculate atmospheric boundary layer height from the inversion of + absolute temperature gradient + + It is the height where the temperature gradient (absolute or potential) changes of sign. + See [Col14]. + + Parameters + ---------- + height : `pint.Quantity` + Altitude (metres above ground) of the points in the profile + humidity : `pint.Quantity` + Temperature (absolute or potential) profile + smoothingspan : int, optional + The amount of smoothing (number of points in moving average) + idxfoot : int, optional + The index of the foot point (first trusted measure), defaults to 0. + + Returns + ------- + blh : `pint.Quantity` + Boundary layer height estimation + """ + + dTdz = mpcalc.first_derivative(smooth(temperature, smoothingspan), x=height) + + dTdz = dTdz[idxfoot:] + height = height[idxfoot:] + + if any(dTdz * dTdz[0] < 0): + iblh = np.where(dTdz * dTdz[0] < 0)[0][0] + blh = height[iblh] + else: + blh = np.nan * units.meter + + return blh diff --git a/tests/test_boundarylayer.py b/tests/test_boundarylayer.py new file mode 100644 index 0000000000..d09e60af32 --- /dev/null +++ b/tests/test_boundarylayer.py @@ -0,0 +1,145 @@ +#!/usr/bin/python +# -*-coding:utf-8 -*- +"""Testing program for the MetPy boundary layer module""" + +import numpy as np +from metpy.units import units +from metpy.calc import boundarylayer +import metpy.calc as mpcalc + +# SAMPLE DATA +# =========== + +# Raw data +# -------- +# Sample taken from a real radiosonding made at Trappes (France) the 2018-08-02 11:01 + +SAMPLE_HEIGHT = np.array( + [ + 189.95, 228.08, 282.48, 337.66, 380.98, 424.59, 456.27, 490.51, + 525.4, 561.28, 595.6, 620.15, 654.14, 697.99, 742.43, 778.43, + 820.95, 870.18, 916.58, 960.31, 1004.73, 1050.9, 1109.05, 1162.18, + 1203.71, 1255.9, 1302.08, 1345.05, 1393.14, 1444.82, 1485.64, 1541.92, + 1590.16, 1631.41, 1674.2, 1716.32, 1767.26, 1818.96, 1858.13, 1905.87, + 1955.06, 2005.74, 2051.76, 2102.09, 2154.78, 2201.7, 2254.76, 2300.89, + 2345.48, 2404.11, 2448.66, 2502.18, 2561.29, 2603.77, 2652.56, + 2704.48, 2747.04, 2794.74, 2841.13, 2891.76, 2938.83, 2991.56, + 3041.57, 3104.81, 3161.18, 3219.26, 3268.16, 3312.8, 3365.91, 3403.13, + 3451.95, 3511.36, 3561.1, 3610.11, 3647.65, 3687.3, 3735.96, 3801.44, + 3849.53, 3908.12, 3965.21, 4016.0, 4062.18, 4104.46, 4153.43, 4191.55, + 4237.61, 4280.69, 4317.76, 4359.39, 4405.79, 4445.88, 4490.97, 4539.87, + 4584.26, 4635.39, 4677.75, 4714.34, 4753.45, 4793.44 + ] +) * units.metres + +SAMPLE_TEMPERATURE = np.array( + [ + 298.21, 299.87, 299.16, 298.61, 298.15, 297.7 , 297.37, 297.01, + 296.57, 296.25, 295.93, 295.64, 295.36, 294.99, 294.54, 294.17, + 293.78, 293.29, 292.9 , 292.41, 292.01, 291.55, 291.04, 290.8 , + 290.48, 290.04, 289.67, 289.26, 288.77, 288.33, 287.95, 287.42, + 286.94, 286.54, 286.17, 285.8 , 285.36, 284.99, 284.72, 284.33, + 283.85, 283.82, 284.01, 284.06, 285.07, 285.63, 285.7 , 285.58, + 285.31, 284.99, 284.91, 284.72, 284.5 , 284.33, 284.16, 283.94, + 283.67, 283.43, 283.19, 282.81, 282.48, 282.2 , 281.96, 281.79, + 281.41, 280.92, 280.59, 280.17, 279.68, 279.38, 279.01, 278.71, + 278.44, 278.43, 278.36, 278.25, 278.16, 277.54, 277.02, 276.45, + 276.4 , 276.42, 276.26, 275.99, 275.65, 275.36, 274.97, 274.7 , + 274.5 , 274.29, 274.07, 273.88, 273.73, 273.46, 273.21, 272.78, + 272.49, 272.33, 272.24, 272.46 + ] +) * units.kelvin + +SAMPLE_PRESSURE = np.array( + [ + 1002.2, 997.9, 991.8, 985.6, 980.7, 975.8, 972.3, 968.5, + 964.7, 960.7, 956.9, 954.2, 950.5, 945.7, 940.9, 937. , + 932.4, 927.1, 922.1, 917.5, 912.7, 907.8, 901.7, 896.1, + 891.8, 886.4, 881.6, 877.1, 872.2, 866.9, 862.7, 857. , + 852.1, 848. , 843.7, 839.5, 834.4, 829.3, 825.4, 820.7, + 815.9, 811. , 806.5, 801.6, 796.6, 792.1, 787.1, 782.8, + 778.7, 773.2, 769.1, 764.2, 758.8, 754.9, 750.5, 745.9, + 742.1, 737.8, 733.7, 729.3, 725.1, 720.5, 716.2, 710.7, + 705.9, 700.9, 696.8, 693. , 688.6, 685.4, 681.4, 676.4, + 672.3, 668.3, 665.2, 662. , 658.1, 652.8, 649. , 644.3, + 639.8, 635.8, 632.1, 628.8, 625. , 622.1, 618.6, 615.3, + 612.4, 609.3, 605.7, 602.7, 599.3, 595.7, 592.4, 588.6, + 585.5, 582.8, 580. , 577.1 + ] +) * units.hPa + +SAMPLE_U = np.array( + [ + 0. , -0.73, -1.5 , -0.59, -4.14, -4.51, -4.43, 2.96, -3.86, + -2.88, -5.79, -4.24, 0.58, -2.2 , 1.32, -7.32, -3.77, 2.55, + 0.72, -1.2 , -6.18, 1.37, 1.58, 0.8 , -3.12, -0.15, -4.4 , + -1.47, -4.68, 1.28, 5.19, -3. , -4.52, -2.9 , 3.09, 4.48, + -0.15, -2.97, -0.66, 0.57, 2.18, 0.31, -3.66, 1.06, 5.57, + -3.97, 1.53, 1.78, -3.44, -2.01, -0.88, 1.65, -0.61, 2.06, + 4.27, 0.69, -3.39, 2.82, -0.58, -3.63, -2.08, -4.77, -4.05, + -3.07, 3.34, -6.44, -0.3 , 1.55, -3.61, -0.01, 4. , -4.53, + 5.44, 3.21, 2.98, 4.87, 2.02, 2.42, 3.78, 0.89, -5.65, + -1.22, 4.81, -0.24, -0.72, 1.79, 3. , 2.38, 1.38, 6.35, + 9.45, 3.71, -1.24, 2.25, 2.12, 2.8 , 5.71, 5.61, 6.78, + 7.99 + ] +) * units.meter_per_second + +SAMPLE_V = np.array( + [ + -0.01, -1.56, -3.73, 7.87, -2.07, -1.8, -5.91, 3.86, -1.62, -2.82, + -2.19, -4.74, -5.87, 0.32, -3.34, -0.71, 4.32, -0.97, 3.56, 5.26, + 2.46, -1.88, -4.58, -3.45, 2.76, -3.63, 1.13, 4.07, -2.67, -4.55, + 2.43, -0.38, -1.97, -0.58, -2.32, -1.24, 2.76, -1.25, -2.11, 3.64, + -3.37, -3.37, 3.03, 4.87, -0.88, 0.25, 0.52, -6.89, -3.69, -0.01, + 1.7, 0.47, 0.46, 0.03, 0.17, 1.4, 1.97, -1.81, 0.1, -5.69, -5.77, + -6.36, -6.09, -0.36, -1.96, -5.48, 1.39, -0.93, -3.39, -0.04, 1.71, + 1.68, 1.59, -4.38, -2.79, -4.65, -1.84, -0.01, 3.58, -3.44, -4.42, + -8.27, -11.55, -10.06, -9.23, -9.33, -4.06, -5.7, -5.42, -5.49, + -6.71, -7.8, -4.16, -3.03, -8.58, -7.86, -4.35, -4.4, -4.45, -1.89 + ] +) * units.meter_per_second + +SAMPLE_RELATIVE_HUMIDITY = np.array( + [ + 48.5, 45.2, 45. , 46. , 47.6, 48.2, 49.4, 50.5, 51.6, 52.4, 53.2, + 54.4, 55.3, 55.5, 56.7, 58.2, 59. , 60.6, 61.2, 62.9, 64.4, 65.7, + 66.7, 67.1, 66.6, 66.9, 67. , 67.7, 69.6, 71.6, 72.8, 73.3, 75.5, + 76.9, 78.1, 79.4, 81.3, 82.1, 80.6, 80.8, 82.9, 76.9, 62.1, 57.6, + 29.1, 24.4, 21.9, 28.3, 28.5, 28.2, 28.1, 27.1, 28.7, 28.7, 25.8, + 31.3, 31.6, 33.9, 32.4, 37. , 40. , 40.9, 37.7, 37.2, 41. , 43.4, + 44.8, 49. , 50.6, 52.8, 49.9, 46.4, 39.4, 28.3, 30.7, 26.7, 23.8, + 41.3, 45. , 48.2, 27.7, 18. , 14.8, 16.1, 17.3, 16.4, 16.2, 15.7, + 14.9, 14.8, 14.5, 13.4, 12.4, 14.6, 15.3, 15.5, 15.3, 14.7, 13.8, + 11.7 + ] +) * units.percent + + +# Deduced data +# ------------ +potential_temperature = mpcalc.potential_temperature(SAMPLE_PRESSURE, SAMPLE_TEMPERATURE) +dewpoint = mpcalc.dewpoint_from_relative_humidity(SAMPLE_TEMPERATURE, SAMPLE_RELATIVE_HUMIDITY) +specific_humidity = mpcalc.specific_humidity_from_dewpoint(SAMPLE_PRESSURE, dewpoint) + + +# BOUNDARY LAYER HEIGHT ESTIMATIONS +# ================================= + +blh = boundarylayer.blh_from_richardson_bulk(SAMPLE_HEIGHT, potential_temperature, SAMPLE_U, SAMPLE_V) +print(f"Estimated boundary layer height with blh_from_richardson_bulk:\t {blh}") + +blh = boundarylayer.blh_from_parcel(SAMPLE_HEIGHT, potential_temperature) +print(f"Estimated boundary layer height with blh_from_parcel:\t {blh}") + +blh = boundarylayer.blh_from_humidity_gradient(SAMPLE_HEIGHT, SAMPLE_RELATIVE_HUMIDITY) +print(f"Estimated boundary layer height with blh_from_humidity_gradient (relative humidity):\t {blh}") + +blh = boundarylayer.blh_from_humidity_gradient(SAMPLE_HEIGHT, specific_humidity) +print(f"Estimated boundary layer height with blh_from_humidity_gradient (specific humidity):\t {blh}") + +blh = boundarylayer.blh_from_temperature_inversion(SAMPLE_HEIGHT, SAMPLE_TEMPERATURE) +print(f"Estimated boundary layer height with blh_from_temperature_inversion (absolute temperature):\t {blh}") + +blh = boundarylayer.blh_from_temperature_inversion(SAMPLE_HEIGHT, potential_temperature) +print(f"Estimated boundary layer height with blh_from_temperature_inversion (potential temperature):\t {blh}") From 6bdf81bae326153a465d3ac1485b0628e29ec456 Mon Sep 17 00:00:00 2001 From: Thomas Rieutord Date: Wed, 31 Jul 2024 13:03:38 +0100 Subject: [PATCH 2/5] Update bulk_richardson_number: fix null wind at bottom (comment #2 in PR #3572) --- src/metpy/calc/boundarylayer.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/metpy/calc/boundarylayer.py b/src/metpy/calc/boundarylayer.py index 1d43ea0ebe..673b3f893e 100644 --- a/src/metpy/calc/boundarylayer.py +++ b/src/metpy/calc/boundarylayer.py @@ -93,13 +93,15 @@ def bulk_richardson_number( `pint.Quantity` Bulk Richardson number profile """ - - u[0] = 0 * units.meter_per_second - v[0] = 0 * units.meter_per_second - + if idxfoot == 0: + # Force the ground level to have null wind + Du = u + Dv = v + else: + Du = u - u[idxfoot] + Dv = v - v[idxfoot] + Dtheta = potential_temperature - potential_temperature[idxfoot] - Du = u - u[idxfoot] - Dv = v - v[idxfoot] Dz = height - height[idxfoot] idx0 = Du**2 + Dv**2 + ustar**2 == 0 From 106ac48ccd5e7a63aaa55714b5f328702311fec7 Mon Sep 17 00:00:00 2001 From: Thomas Rieutord Date: Wed, 31 Jul 2024 13:14:19 +0100 Subject: [PATCH 3/5] Rename blh_from_humidity_gradient -> blh_from_concentration_gradient --- src/metpy/calc/boundarylayer.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/metpy/calc/boundarylayer.py b/src/metpy/calc/boundarylayer.py index 673b3f893e..8ed9d68a81 100644 --- a/src/metpy/calc/boundarylayer.py +++ b/src/metpy/calc/boundarylayer.py @@ -225,24 +225,24 @@ def blh_from_parcel( return blh -def blh_from_humidity_gradient( +def blh_from_concentration_gradient( height, - humidity, + concentration_profile, smoothingspan: int = 5, idxfoot: int = 0, ): - """Calculate atmospheric boundary layer height from the relative - humidity gradient + """Calculate atmospheric boundary layer height from a concentration + profile (specific/relative humidity, aerosol backscatter, TKE..) - It is the height where the relative humidity or specific humidity gradient reaches a minimum. + It is the height where the gradient of the concentration profile reaches a minimum. See [Sei00, HL06, Col14]. Parameters ---------- height : `pint.Quantity` Altitude (metres above ground) of the points in the profile - humidity : `pint.Quantity` - Humidity (relative or specific) profile + concentration_profile : `pint.Quantity` + Concentration profile (specific/relative humidity, aerosol backscatter, TKE..) smoothingspan : int, optional The amount of smoothing (number of points in moving average) idxfoot : int, optional @@ -254,12 +254,10 @@ def blh_from_humidity_gradient( Boundary layer height estimation """ - dRHdz = mpcalc.first_derivative(smooth(humidity, smoothingspan), x=height) - - dRHdz = dRHdz[idxfoot:] + dcdz = mpcalc.first_derivative(smooth(concentration_profile, smoothingspan), x=height) + dcdz = dcdz[idxfoot:] height = height[idxfoot:] - - iblh = np.argmin(dRHdz) + iblh = np.argmin(dcdz) return height[iblh] From 0f16c2d395a604b27313ce3014c79fd503566565 Mon Sep 17 00:00:00 2001 From: Thomas Rieutord Date: Wed, 31 Jul 2024 14:10:19 +0100 Subject: [PATCH 4/5] Improve doc + linting (comment #3 in PR #3572) --- src/metpy/calc/boundarylayer.py | 48 +++++++++++++++++---------------- 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/src/metpy/calc/boundarylayer.py b/src/metpy/calc/boundarylayer.py index 8ed9d68a81..85e4b0c662 100644 --- a/src/metpy/calc/boundarylayer.py +++ b/src/metpy/calc/boundarylayer.py @@ -1,13 +1,12 @@ -#!/usr/bin/python -# -*-coding:utf-8 -*- +# Copyright (c) 2024 MetPy Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause """ Contains a collection of boundary layer height estimations. - -References: ------------ - -[Col14]: Collaud Coen, M., Praz, C., Haefele, A., Ruffieux, D., Kaufmann, P., and Calpini, B. (2014): +References +---------- +[Col14]: Collaud Coen, M., Praz, C., Haefele, A., Ruffieux, D., Kaufmann, P., and Calpini, B. (2014) Determination and climatology of the planetary boundary layer height above the Swiss plateau by in situ and remote sensing measurements as well as by the COSMO-2 model Atmos. Chem. Phys., 14, 13205–13221. @@ -15,19 +14,18 @@ Determination of the atmospheric boundary layer height from radiosonde and lidar backscatter. Boundary-Layer Meteorology, 120(1), 181-200. -[Guo16]: Guo, J., Miao, Y., Zhang, Y., Liu, H., Li, Z., Zhang, W., ... & Zhai, P. (2016): +[Guo16]: Guo, J., Miao, Y., Zhang, Y., Liu, H., Li, Z., Zhang, W., ... & Zhai, P. (2016) The climatology of planetary boundary layer height in China derived from radiosonde and reanalysis data. Atmos. Chem. Phys, 16(20), 13309-13319. -[Sei00]: Seidel, D. J., Ao, C. O., & Li, K. (2010): +[Sei00]: Seidel, D. J., Ao, C. O., & Li, K. (2010) Estimating climatological planetary boundary layer heights from radiosonde observations: Comparison of methods and uncertainty analysis. Journal of Geophysical Research: Atmospheres, 115(D16). -[VH96]: Vogelezang, D. H. P., & Holtslag, A. A. M. (1996): +[VH96]: Vogelezang, D. H. P., & Holtslag, A. A. M. (1996) Evaluation and model impacts of alternative boundary-layer height formulations. Boundary-Layer Meteorology, 81(3-4), 245-269. """ - import numpy as np from copy import deepcopy @@ -51,11 +49,17 @@ def smooth(val, span): ------- smoothed_val: array-like Array of smoothed values + + See also + -------- + [`bottleneck.move_mean`](https://bottleneck.readthedocs.io/en/latest/reference.html#bottleneck.move_mean), + [`scipy.ndimage.uniform_filter1d`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.uniform_filter1d.html#scipy.ndimage.uniform_filter1d), + [`pandas.DataFrame.rolling`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.rolling.html) """ - N = len(val) + n = len(val) smoothed_val = deepcopy(val) - for i in range(N): - smoothed_val[i] = np.nanmean(val[i - min(span, i) : i + min(span, N - i)]) + for i in range(n): + smoothed_val[i] = np.nanmean(val[i - min(span, i) : i + min(span, n - i)]) return smoothed_val @@ -136,7 +140,7 @@ def blh_from_richardson_bulk( bulk Richardson number. It is the height where the bulk Richardson number exceeds a given threshold. - See [VH96, Sei00, Col14, Guo16]. + Well indicated for unstable boundary layers. See [VH96, Sei00, Col14, Guo16]. Parameters ---------- @@ -189,11 +193,11 @@ def blh_from_parcel( smoothingspan: int = 5, theta0=None, ): - """Calculate atmospheric boundary layer height with the parcel method. + """Calculate atmospheric boundary layer height with the "parcel method" + (or "potential temperature threshold method"). - It is the height where the potential temperature profile reaches its - foot value. - See [Sei00, HL06, Col14]. + It is the height where the potential temperature profile exceeds its + foot value. Well indicated for unstable boundary layers. See [Sei00, HL06, Col14]. Parameters ---------- @@ -235,7 +239,7 @@ def blh_from_concentration_gradient( profile (specific/relative humidity, aerosol backscatter, TKE..) It is the height where the gradient of the concentration profile reaches a minimum. - See [Sei00, HL06, Col14]. + Well indicated for stable boundary layers. See [Sei00, HL06, Col14]. Parameters ---------- @@ -253,7 +257,6 @@ def blh_from_concentration_gradient( blh : `pint.Quantity` Boundary layer height estimation """ - dcdz = mpcalc.first_derivative(smooth(concentration_profile, smoothingspan), x=height) dcdz = dcdz[idxfoot:] height = height[idxfoot:] @@ -272,7 +275,7 @@ def blh_from_temperature_inversion( absolute temperature gradient It is the height where the temperature gradient (absolute or potential) changes of sign. - See [Col14]. + Well indicated for stable boundary layers. See [Col14]. Parameters ---------- @@ -290,7 +293,6 @@ def blh_from_temperature_inversion( blh : `pint.Quantity` Boundary layer height estimation """ - dTdz = mpcalc.first_derivative(smooth(temperature, smoothingspan), x=height) dTdz = dTdz[idxfoot:] From c60a30b623f9c5933c70c24f24ba11608078c4e8 Mon Sep 17 00:00:00 2001 From: Thomas Rieutord Date: Wed, 31 Jul 2024 14:12:10 +0100 Subject: [PATCH 5/5] Rewrite tests with standard data and for pytest (comment #7 in PR #3572) --- tests/test_boundarylayer.py | 171 ++++++++++-------------------------- 1 file changed, 46 insertions(+), 125 deletions(-) diff --git a/tests/test_boundarylayer.py b/tests/test_boundarylayer.py index d09e60af32..98ab7718f9 100644 --- a/tests/test_boundarylayer.py +++ b/tests/test_boundarylayer.py @@ -3,143 +3,64 @@ """Testing program for the MetPy boundary layer module""" import numpy as np -from metpy.units import units -from metpy.calc import boundarylayer +import pandas as pd + import metpy.calc as mpcalc +from metpy.calc import boundarylayer +from metpy.cbook import get_test_data +from metpy.units import units # SAMPLE DATA # =========== +col_names = ["pressure", "height", "temperature", "dewpoint", "direction", "speed"] -# Raw data -# -------- -# Sample taken from a real radiosonding made at Trappes (France) the 2018-08-02 11:01 - -SAMPLE_HEIGHT = np.array( - [ - 189.95, 228.08, 282.48, 337.66, 380.98, 424.59, 456.27, 490.51, - 525.4, 561.28, 595.6, 620.15, 654.14, 697.99, 742.43, 778.43, - 820.95, 870.18, 916.58, 960.31, 1004.73, 1050.9, 1109.05, 1162.18, - 1203.71, 1255.9, 1302.08, 1345.05, 1393.14, 1444.82, 1485.64, 1541.92, - 1590.16, 1631.41, 1674.2, 1716.32, 1767.26, 1818.96, 1858.13, 1905.87, - 1955.06, 2005.74, 2051.76, 2102.09, 2154.78, 2201.7, 2254.76, 2300.89, - 2345.48, 2404.11, 2448.66, 2502.18, 2561.29, 2603.77, 2652.56, - 2704.48, 2747.04, 2794.74, 2841.13, 2891.76, 2938.83, 2991.56, - 3041.57, 3104.81, 3161.18, 3219.26, 3268.16, 3312.8, 3365.91, 3403.13, - 3451.95, 3511.36, 3561.1, 3610.11, 3647.65, 3687.3, 3735.96, 3801.44, - 3849.53, 3908.12, 3965.21, 4016.0, 4062.18, 4104.46, 4153.43, 4191.55, - 4237.61, 4280.69, 4317.76, 4359.39, 4405.79, 4445.88, 4490.97, 4539.87, - 4584.26, 4635.39, 4677.75, 4714.34, 4753.45, 4793.44 - ] -) * units.metres - -SAMPLE_TEMPERATURE = np.array( - [ - 298.21, 299.87, 299.16, 298.61, 298.15, 297.7 , 297.37, 297.01, - 296.57, 296.25, 295.93, 295.64, 295.36, 294.99, 294.54, 294.17, - 293.78, 293.29, 292.9 , 292.41, 292.01, 291.55, 291.04, 290.8 , - 290.48, 290.04, 289.67, 289.26, 288.77, 288.33, 287.95, 287.42, - 286.94, 286.54, 286.17, 285.8 , 285.36, 284.99, 284.72, 284.33, - 283.85, 283.82, 284.01, 284.06, 285.07, 285.63, 285.7 , 285.58, - 285.31, 284.99, 284.91, 284.72, 284.5 , 284.33, 284.16, 283.94, - 283.67, 283.43, 283.19, 282.81, 282.48, 282.2 , 281.96, 281.79, - 281.41, 280.92, 280.59, 280.17, 279.68, 279.38, 279.01, 278.71, - 278.44, 278.43, 278.36, 278.25, 278.16, 277.54, 277.02, 276.45, - 276.4 , 276.42, 276.26, 275.99, 275.65, 275.36, 274.97, 274.7 , - 274.5 , 274.29, 274.07, 273.88, 273.73, 273.46, 273.21, 272.78, - 272.49, 272.33, 272.24, 272.46 - ] -) * units.kelvin - -SAMPLE_PRESSURE = np.array( - [ - 1002.2, 997.9, 991.8, 985.6, 980.7, 975.8, 972.3, 968.5, - 964.7, 960.7, 956.9, 954.2, 950.5, 945.7, 940.9, 937. , - 932.4, 927.1, 922.1, 917.5, 912.7, 907.8, 901.7, 896.1, - 891.8, 886.4, 881.6, 877.1, 872.2, 866.9, 862.7, 857. , - 852.1, 848. , 843.7, 839.5, 834.4, 829.3, 825.4, 820.7, - 815.9, 811. , 806.5, 801.6, 796.6, 792.1, 787.1, 782.8, - 778.7, 773.2, 769.1, 764.2, 758.8, 754.9, 750.5, 745.9, - 742.1, 737.8, 733.7, 729.3, 725.1, 720.5, 716.2, 710.7, - 705.9, 700.9, 696.8, 693. , 688.6, 685.4, 681.4, 676.4, - 672.3, 668.3, 665.2, 662. , 658.1, 652.8, 649. , 644.3, - 639.8, 635.8, 632.1, 628.8, 625. , 622.1, 618.6, 615.3, - 612.4, 609.3, 605.7, 602.7, 599.3, 595.7, 592.4, 588.6, - 585.5, 582.8, 580. , 577.1 - ] -) * units.hPa - -SAMPLE_U = np.array( - [ - 0. , -0.73, -1.5 , -0.59, -4.14, -4.51, -4.43, 2.96, -3.86, - -2.88, -5.79, -4.24, 0.58, -2.2 , 1.32, -7.32, -3.77, 2.55, - 0.72, -1.2 , -6.18, 1.37, 1.58, 0.8 , -3.12, -0.15, -4.4 , - -1.47, -4.68, 1.28, 5.19, -3. , -4.52, -2.9 , 3.09, 4.48, - -0.15, -2.97, -0.66, 0.57, 2.18, 0.31, -3.66, 1.06, 5.57, - -3.97, 1.53, 1.78, -3.44, -2.01, -0.88, 1.65, -0.61, 2.06, - 4.27, 0.69, -3.39, 2.82, -0.58, -3.63, -2.08, -4.77, -4.05, - -3.07, 3.34, -6.44, -0.3 , 1.55, -3.61, -0.01, 4. , -4.53, - 5.44, 3.21, 2.98, 4.87, 2.02, 2.42, 3.78, 0.89, -5.65, - -1.22, 4.81, -0.24, -0.72, 1.79, 3. , 2.38, 1.38, 6.35, - 9.45, 3.71, -1.24, 2.25, 2.12, 2.8 , 5.71, 5.61, 6.78, - 7.99 - ] -) * units.meter_per_second - -SAMPLE_V = np.array( - [ - -0.01, -1.56, -3.73, 7.87, -2.07, -1.8, -5.91, 3.86, -1.62, -2.82, - -2.19, -4.74, -5.87, 0.32, -3.34, -0.71, 4.32, -0.97, 3.56, 5.26, - 2.46, -1.88, -4.58, -3.45, 2.76, -3.63, 1.13, 4.07, -2.67, -4.55, - 2.43, -0.38, -1.97, -0.58, -2.32, -1.24, 2.76, -1.25, -2.11, 3.64, - -3.37, -3.37, 3.03, 4.87, -0.88, 0.25, 0.52, -6.89, -3.69, -0.01, - 1.7, 0.47, 0.46, 0.03, 0.17, 1.4, 1.97, -1.81, 0.1, -5.69, -5.77, - -6.36, -6.09, -0.36, -1.96, -5.48, 1.39, -0.93, -3.39, -0.04, 1.71, - 1.68, 1.59, -4.38, -2.79, -4.65, -1.84, -0.01, 3.58, -3.44, -4.42, - -8.27, -11.55, -10.06, -9.23, -9.33, -4.06, -5.7, -5.42, -5.49, - -6.71, -7.8, -4.16, -3.03, -8.58, -7.86, -4.35, -4.4, -4.45, -1.89 - ] -) * units.meter_per_second - -SAMPLE_RELATIVE_HUMIDITY = np.array( - [ - 48.5, 45.2, 45. , 46. , 47.6, 48.2, 49.4, 50.5, 51.6, 52.4, 53.2, - 54.4, 55.3, 55.5, 56.7, 58.2, 59. , 60.6, 61.2, 62.9, 64.4, 65.7, - 66.7, 67.1, 66.6, 66.9, 67. , 67.7, 69.6, 71.6, 72.8, 73.3, 75.5, - 76.9, 78.1, 79.4, 81.3, 82.1, 80.6, 80.8, 82.9, 76.9, 62.1, 57.6, - 29.1, 24.4, 21.9, 28.3, 28.5, 28.2, 28.1, 27.1, 28.7, 28.7, 25.8, - 31.3, 31.6, 33.9, 32.4, 37. , 40. , 40.9, 37.7, 37.2, 41. , 43.4, - 44.8, 49. , 50.6, 52.8, 49.9, 46.4, 39.4, 28.3, 30.7, 26.7, 23.8, - 41.3, 45. , 48.2, 27.7, 18. , 14.8, 16.1, 17.3, 16.4, 16.2, 15.7, - 14.9, 14.8, 14.5, 13.4, 12.4, 14.6, 15.3, 15.5, 15.3, 14.7, 13.8, - 11.7 - ] -) * units.percent - - -# Deduced data -# ------------ -potential_temperature = mpcalc.potential_temperature(SAMPLE_PRESSURE, SAMPLE_TEMPERATURE) -dewpoint = mpcalc.dewpoint_from_relative_humidity(SAMPLE_TEMPERATURE, SAMPLE_RELATIVE_HUMIDITY) -specific_humidity = mpcalc.specific_humidity_from_dewpoint(SAMPLE_PRESSURE, dewpoint) +df = pd.read_fwf( + get_test_data("may4_sounding.txt", as_file_obj=False), + skiprows=5, + usecols=[0, 1, 2, 3, 6, 7], + names=col_names, +) + +# Drop any rows with all NaN values for T, Td, winds +df = df.dropna( + subset=("temperature", "dewpoint", "direction", "speed"), how="all" +).reset_index(drop=True) + +height = df["height"].values * units.metres +pressure = df["pressure"].values * units.hPa +temperature = df["temperature"].values * units.degC +dewpoint = df["dewpoint"].values * units.degC +wind_speed = df["speed"].values * units.knots +wind_dir = df["direction"].values * units.degrees + +u, v = mpcalc.wind_components(wind_speed, wind_dir) +relative_humidity = mpcalc.relative_humidity_from_dewpoint(temperature, dewpoint) +potential_temperature = mpcalc.potential_temperature(pressure, temperature) +specific_humidity = mpcalc.specific_humidity_from_dewpoint(pressure, dewpoint) # BOUNDARY LAYER HEIGHT ESTIMATIONS # ================================= -blh = boundarylayer.blh_from_richardson_bulk(SAMPLE_HEIGHT, potential_temperature, SAMPLE_U, SAMPLE_V) -print(f"Estimated boundary layer height with blh_from_richardson_bulk:\t {blh}") +def test_blh_from_richardson_bulk(): + blh = boundarylayer.blh_from_richardson_bulk(height, potential_temperature, u, v) + blh_true = 1397 * units.meter + assert blh == blh_true + -blh = boundarylayer.blh_from_parcel(SAMPLE_HEIGHT, potential_temperature) -print(f"Estimated boundary layer height with blh_from_parcel:\t {blh}") +def test_blh_from_parcel(): + blh = boundarylayer.blh_from_parcel(height, potential_temperature) + blh_true = 610 * units.meter + assert blh == blh_true -blh = boundarylayer.blh_from_humidity_gradient(SAMPLE_HEIGHT, SAMPLE_RELATIVE_HUMIDITY) -print(f"Estimated boundary layer height with blh_from_humidity_gradient (relative humidity):\t {blh}") -blh = boundarylayer.blh_from_humidity_gradient(SAMPLE_HEIGHT, specific_humidity) -print(f"Estimated boundary layer height with blh_from_humidity_gradient (specific humidity):\t {blh}") +def test_blh_from_concentration_gradient(): + blh = boundarylayer.blh_from_concentration_gradient(height, specific_humidity) + blh_true = 1766 * units.meter + assert blh == blh_true -blh = boundarylayer.blh_from_temperature_inversion(SAMPLE_HEIGHT, SAMPLE_TEMPERATURE) -print(f"Estimated boundary layer height with blh_from_temperature_inversion (absolute temperature):\t {blh}") -blh = boundarylayer.blh_from_temperature_inversion(SAMPLE_HEIGHT, potential_temperature) -print(f"Estimated boundary layer height with blh_from_temperature_inversion (potential temperature):\t {blh}") +def test_blh_from_temperature_inversion(): + blh = boundarylayer.blh_from_temperature_inversion(height, potential_temperature) + blh_true = 610 * units.meter + assert blh == blh_true