Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.2k stars 405 forks source link

Support for Gridded data #3481

Open leaver2000 opened 1 month ago

leaver2000 commented 1 month ago

What should we add?

I've been developing an application that works with gridded data. The metpy.calc.thermo module has been a great guide for my work thus far.

My moist_lapse function differs from metpy.calc.thermo.moist_lapse and is written in Cython, that seemed to be one of the major bottlenecks the vectorizing things for gridded data support.

I currently have implementations of the moist_lapse, wet_bulb_temperature, parcel_profile, ccl, and downdraft_cape that support 2d: (N, Z) array structure.

My implementations don't use pint, everything is assumed si units. Let me know if there is any interest I'm happy to share.

pressure = isobaric_levels()
(temperature, specific_humidity), shape = get_values(["temperatureair", "humidityspecific"])
dewpoint = thermo.dewpoint_from_specific_humidity(pressure[newaxis, :], specific_humidity)

dcape = thermo.downdraft_cape(pressure, temperature, dewpoint)
dcape = dcape.reshape((-1,) + shape[2:])

fig, axes = plt.subplots(1, dcape.shape[0])
fig.tight_layout()
fig.suptitle("dcape", fontsize=16, x=0.5, y=0.65)

for i, ax in enumerate(axes):
    ax.imshow(dcape[i, :, :], cmap="coolwarm")

image

T = temperature[::2500]
Td = dewpoint[::2500]
my_dcape = thermo.downdraft_cape(
    pressure=pressure,
    temperature=T,
    dewpoint=Td,
)
mp_dcape = np.array(
    [
        mpcalc.downdraft_cape(pressure * units.pascal, T[i] * units.kelvin, Td[i] * units.kelvin)[0].m
        for i in range(T.shape[0])
    ]
)

delta = np.abs(my_dcape - mp_dcape)
print(f"{delta.max() = } {delta.mean() = } {delta.std() = }")

delta.max() = 0.6819282354549614 delta.mean() = 0.1615940182263646 delta.std() = 0.13604876693905152

Reference

No response

leaver2000 commented 1 month ago
assert (260640, 24) == temperature.shape == dewpoint.shape # (T*Y*X, Z)
assert (24,) == pressure.shape
%timeit thermo.downdraft_cape(pressure, temperature, dewpoint)
# 218 ms ± 11 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Z-Richard commented 1 month ago

First of all, great work! Indeed, there has been considerable debate and effort on making MetPy faster through (a) vectorization, e.g., #1968; and (b) using faster variants of NumPy, e.g., #3432. I think this increasing need echoes the demand of calculating convective metrics (CAPE, CIN) with model data, as understanding the changing nature of thunderstorms in a warming world has become a very important research topic. The original development of many MetPy thermodynamic functions, however, likely has a meteorological audience in mind (who might only need to calculate CAPE/CIN at a few points). And yet, the scientific community currently lacks such a tool to calculate these quantities en masse. It would be helpful (and deeply appreciated) if you can put your code in a public repository (as a reference point if the MetPy developers want to take a look at the code later) or directly submit a draft PR (like #1968). While I am not a core developer of MetPy and do not know how this effort will evolve, I believe that building a faster, vectorized version of thermodynamic functions as a community will benefit future users substantially.

leaver2000 commented 1 month ago

@dopplershift what type of precision are you looking for gridded solutions? I've put together a couple of benchmarks below.

Cython functions

moist_lapse

import numpy as np
import metpy.calc as mpcalc
from metpy.units import units

import nzthermo as nzt
N = 1000
Z = 20

P = np.linspace(101325, 10000, Z)[np.newaxis, :] # (1, Z)
T = np.random.uniform(300, 200, N) # (N,)

ml = nzt.moist_lapse(P, T)
%timeit nzt.moist_lapse(P, T)
# 1.22 ms ± 142 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
P = P[0] * units.Pa
T = T * units.kelvin
ml_ = [mpcalc.moist_lapse(P, T[i]).m for i in range(N)]  # type: ignore
%timeit [mpcalc.moist_lapse(P, T[i]).m for i in range(N)]
# 1.65 s ± 29.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
np.testing.assert_allclose(ml, ml_, rtol=1e-3)

lcl

P = np.random.uniform(101325, 10000, N) # (N,)
T = np.random.uniform(300, 200, N) # (N,)
Td = T - np.random.uniform(0, 10, N) # (N,)

lcl_p, lcl_t = nzt.lcl(P, T, Td) # ((N,), (N,))
%timeit nzt.lcl(P, T, Td)
# 1.4 ms ± 373 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
P *= units.Pa
T *= units.kelvin
Td *= units.kelvin
lcl_p_, lcl_t_ = (x.m for x in mpcalc.lcl(P, T, Td))  # type: ignore
%timeit mpcalc.lcl(P, T, Td)
# 1.57 s ± 7.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
np.testing.assert_allclose(lcl_p, lcl_p_, rtol=1e-3)
np.testing.assert_allclose(lcl_t, lcl_t_, rtol=1e-3)

Implementations

wet_bulb_temperature

P = np.random.uniform(101325, 10000, 1000).astype(np.float32)
T = np.random.uniform(300, 200, 1000).astype(np.float32)
Td = T - np.random.uniform(0, 10, 1000).astype(np.float32)

wb = nzt.wet_bulb_temperature(P, T, Td)
%timeit nzt.wet_bulb_temperature(P, T, Td)
# 1.17 ms ± 124 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
P *= units.Pa
T *= units.kelvin
Td *= units.kelvin
wb_ = mpcalc.wet_bulb_temperature(P, T, Td).m
%timeit mpcalc.wet_bulb_temperature(P, T, Td)
# 390 ms ± 2.86 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
np.testing.assert_allclose(wb, wb_, rtol=1e-3)

downdraft_cape

NOTE: random values caused metpy's dcape function to throw interpolation warnings and return nan in many most cases.

pressure = np.array(
    [1013, 1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 725, 700, 650, 600, 550, 500, 450, 400, 350, 300],
).astype(np.float32)
pressure *= 100.0
temperature = np.array(
    [
        [243, 242, 241, 240, 239, 237, 236, 235, 233, 232, 231, 229, 228, 226, 235, 236, 234, 231, 226, 221, 217, 211],
        [250, 249, 248, 247, 246, 244, 243, 242, 240, 239, 238, 236, 235, 233, 240, 239, 236, 232, 227, 223, 217, 211],
        [293, 292, 290, 288, 287, 285, 284, 282, 281, 279, 279, 280, 279, 278, 275, 270, 268, 264, 260, 254, 246, 237],
        [300, 299, 297, 295, 293, 291, 292, 291, 291, 289, 288, 286, 285, 285, 281, 278, 273, 268, 264, 258, 251, 242],
    ],
    dtype=np.float32,
)
dewpoint = np.array(
    [
        [224, 224, 224, 224, 224, 223, 223, 223, 223, 222, 222, 222, 221, 221, 233, 233, 231, 228, 223, 218, 213, 207],
        [233, 233, 232, 232, 232, 232, 231, 231, 231, 231, 230, 230, 230, 229, 237, 236, 233, 229, 223, 219, 213, 207],
        [288, 288, 287, 286, 281, 280, 279, 277, 276, 275, 270, 258, 244, 247, 243, 254, 262, 248, 229, 232, 229, 224],
        [294, 294, 293, 292, 291, 289, 285, 282, 280, 280, 281, 281, 278, 274, 273, 269, 259, 246, 240, 241, 226, 219],
    ],
    dtype=np.float32,
)

dcape = nzt.downdraft_cape(pressure, temperature, dewpoint)
%timeit nzt.downdraft_cape(pressure, temperature, dewpoint)
# 2.41 ms ± 877 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
P = pressure * units.Pa
T = temperature * units.kelvin
Td = dewpoint * units.kelvin
dcape_ = [mpcalc.downdraft_cape(P, T[i], Td[i])[0].m for i in range(temperature.shape[0])]  # type: ignore
%timeit [mpcalc.downdraft_cape(P, T[i], Td[i]) for i in range(temperature.shape[0])]
# 16.5 ms ± 114 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
np.testing.assert_allclose(dcape, dcape_, rtol=1e-2)

I am currently looking into the _parcel_profile_helper and have a rough implementation developed.

zhangslTHU commented 3 weeks ago

hi, @leaver2000 , thanks a lot for your great work. From above examples, it seems that your implementations are much faster than metpy which takes a long long time to calculate wet bulb temperature.

From the codes and nzthermo website, it seems nzthermo only supports 1D data when calculating wet bulb temperature, could nzthermo calculate 3D data (e.g., time, lat, lon) ?

Last but not least, it would be very helpful if you can provide a way to install nzthermo package with conda (like conda install nzthermo).

leaver2000 commented 3 weeks ago

From the codes and nzthermo website, it seems nzthermo only supports 1D data when calculating wet bulb temperature, could nzthermo calculate 3D data (e.g., time, lat, lon) ?

For wetbulb temperature that can be accomplished as shown below.

TIME = 144
LAT = 721
LON = 1440
N = TIME * LAT * LON
P = np.random.uniform(101325, 10000, N).astype(np.float32).reshape(TIME, LAT, LON)
T = np.random.uniform(300, 200, N).astype(np.float32).reshape(TIME, LAT, LON)
Td = T - np.random.uniform(0, 10, N).astype(np.float32).reshape(TIME, LAT, LON)

wb = nzt.wet_bulb_temperature(P.ravel(), T.ravel(), Td.ravel()).reshape(TIME, LAT, LON)

Last but not least, it would be very helpful if you can provide a way to install nzthermo package with conda (like conda install nzthermo).

I can look into it publishing a release, the challenge is the only supported/tested OS is linux.

zhangslTHU commented 3 weeks ago

Thanks a lot. Hope we can see the released version soon.

How to install nzthermo in linux? I'm curious about the code pip install . shown in nzthermo Getting Started, may be pip install nzthermo ?

leaver2000 commented 3 weeks ago

Thanks a lot. Hope we can see the released version soon.

How to install nzthermo in linux? I'm curious about the code pip install . shown in nzthermo Getting Started, may be pip install nzthermo ?

@zhangslTHU You can use pip to and git to build from source, like...

➜  ~ mkdir myproject
➜  ~ cd myproject 
➜  myproject python3.11 -m venv .venv && source .venv/bin/activate  
(.venv) ➜  myproject pip install git+https://github.com/leaver2000/nzthermo@master
(.venv) ➜  myproject python -c 'import numpy as np;import nzthermo as nzt;print(nzt.wet_bulb_temperature(np.array([101325.]),np.array([300.0]),np.array([292.0])))'
[294.36783585]

pip install . assumes you are inside of the project directory

zhangslTHU commented 3 weeks ago

@leaver2000 thanks a lot. it works now!

It is much faster than Metpy according to test results. It took about 21 mins by Metpy while 0.5s by nzthermo for a array with shape (time=24, lat=181, lon=360),.

Here shows the actual time i spent to calculate wet bulb temperature for ERA5 hourly data in Jan 2022 (31 days) with Metpy, and same data with nzthermo but for 1 year (12 months for 2022):

running time for metpy ``` read 1th data: 2024-05-06 21:36 finish 1th wbt calculation: 2024-05-06 21:57 read 2th data: 2024-05-06 21:57 finish 2th wbt calculation: 2024-05-06 22:18 read 3th data: 2024-05-06 22:18 finish 3th wbt calculation: 2024-05-06 22:40 read 4th data: 2024-05-06 22:40 finish 4th wbt calculation: 2024-05-06 23:01 read 5th data: 2024-05-06 23:01 finish 5th wbt calculation: 2024-05-06 23:22 read 6th data: 2024-05-06 23:22 finish 6th wbt calculation: 2024-05-06 23:44 read 7th data: 2024-05-06 23:44 finish 7th wbt calculation: 2024-05-07 00:05 read 8th data: 2024-05-07 00:05 finish 8th wbt calculation: 2024-05-07 00:27 read 9th data: 2024-05-07 00:27 finish 9th wbt calculation: 2024-05-07 00:48 read 10th data: 2024-05-07 00:48 finish 10th wbt calculation: 2024-05-07 01:09 read 11th data: 2024-05-07 01:09 finish 11th wbt calculation: 2024-05-07 01:31 read 12th data: 2024-05-07 01:31 finish 12th wbt calculation: 2024-05-07 01:52 read 13th data: 2024-05-07 01:52 finish 13th wbt calculation: 2024-05-07 02:14 read 14th data: 2024-05-07 02:14 finish 14th wbt calculation: 2024-05-07 02:35 read 15th data: 2024-05-07 02:35 finish 15th wbt calculation: 2024-05-07 02:56 read 16th data: 2024-05-07 02:56 finish 16th wbt calculation: 2024-05-07 03:18 read 17th data: 2024-05-07 03:18 finish 17th wbt calculation: 2024-05-07 03:39 read 18th data: 2024-05-07 03:39 finish 18th wbt calculation: 2024-05-07 04:00 read 19th data: 2024-05-07 04:00 finish 19th wbt calculation: 2024-05-07 04:21 read 20th data: 2024-05-07 04:21 finish 20th wbt calculation: 2024-05-07 04:43 read 21th data: 2024-05-07 04:43 finish 21th wbt calculation: 2024-05-07 05:04 read 22th data: 2024-05-07 05:04 finish 22th wbt calculation: 2024-05-07 05:25 read 23th data: 2024-05-07 05:25 finish 23th wbt calculation: 2024-05-07 05:47 read 24th data: 2024-05-07 05:47 finish 24th wbt calculation: 2024-05-07 06:08 read 25th data: 2024-05-07 06:08 finish 25th wbt calculation: 2024-05-07 06:29 read 26th data: 2024-05-07 06:29 finish 26th wbt calculation: 2024-05-07 06:50 read 27th data: 2024-05-07 06:50 finish 27th wbt calculation: 2024-05-07 07:12 read 28th data: 2024-05-07 07:12 finish 28th wbt calculation: 2024-05-07 07:33 read 29th data: 2024-05-07 07:33 finish 29th wbt calculation: 2024-05-07 07:54 read 30th data: 2024-05-07 07:54 finish 30th wbt calculation: 2024-05-07 08:16 read 31th data: 2024-05-07 08:16 finish 31th wbt calculation: 2024-05-07 08:37 ````
running time for nzthermo ``` read 1th mon data: 2024-05-07 14:35 finish 1th mon wbt calculation: 2024-05-07 14:35 read 2th mon data: 2024-05-07 14:35 finish 2th mon wbt calculation: 2024-05-07 14:35 read 3th mon data: 2024-05-07 14:35 finish 3th mon wbt calculation: 2024-05-07 14:36 read 4th mon data: 2024-05-07 14:36 finish 4th mon wbt calculation: 2024-05-07 14:36 read 5th mon data: 2024-05-07 14:36 finish 5th mon wbt calculation: 2024-05-07 14:36 read 6th mon data: 2024-05-07 14:36 finish 6th mon wbt calculation: 2024-05-07 14:37 read 7th mon data: 2024-05-07 14:37 finish 7th mon wbt calculation: 2024-05-07 14:37 read 8th mon data: 2024-05-07 14:37 finish 8th mon wbt calculation: 2024-05-07 14:37 read 9th mon data: 2024-05-07 14:37 finish 9th mon wbt calculation: 2024-05-07 14:38 read 10th mon data: 2024-05-07 14:38 finish 10th mon wbt calculation: 2024-05-07 14:38 read 11th mon data: 2024-05-07 14:38 finish 11th mon wbt calculation: 2024-05-07 14:39 read 12th mon data: 2024-05-07 14:39 finish 12th mon wbt calculation: 2024-05-07 14:39 ```

thanks again for the great work and the package your developed, it is very helpful and save much time!

DWesl commented 3 weeks ago

For the functions that require a profile, the signature you chose suggests NumPy gufuncs, have you looked into those? (They are possible in Cython, but dealing with the possibility of someone having a record array of 32-bit floats and 8-bit integers and passing just the floats to your function is annoying).

I'm pretty sure wet-bulb temperature is a point-by-point function (even if you use a definition with parcels rising and falling), so the cython.ufunc decorator will take care of that for you. (Or, since you already have that working, speed up the adaptation of other functions).

leaver2000 commented 3 weeks ago

I've have tested out the cython.ufunc decorator in the past, and I'll admit my current implementations are somewhat inconsistent with typical array broadcasting.

The wetbulb temperature is element-wise ie: (N,)(N,)(N,) but still calls into the moist_lapse ODE which even for a single level is still solved via recursive iteration.

Many other use cases of moist_lapse require (N,Z),(N,),(N,) and it gets a bit tricky when you need each parcel to start and stop at a specific level along Z

dopplershift commented 3 weeks ago

@leaver2000 Apologies for not responding sooner. I want to give a more thoughtful response, but just wanted to say that I am very interested in this work and it aligns with what is our highest development priority, getting things fast enough to have thermo (esp. CAPE) work on input grids.

Regarding precision, I just want to feel confident that what's coming out is correct.

leaver2000 commented 3 weeks ago

@dopplershift cape_cin requires (el, lfc) which are both dependent on the parcel_profile. I've put together a couple implementations of the parcel_profile function. The challenge is because each profile is independent of the next, the return value must be an array of shape (N, Z) where values along Z must be either masked or NaN. This in turn breaks any function downstream of the original implementation. Here is that implementation.

Python 3.11.9 (main, Apr  6 2024, 17:59:24) [GCC 11.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from typing import NamedTuple, Generic, TypeVar, Annotated
>>> import numpy as np
>>> class ParcelProfile(NamedTuple, Generic[float_]):
...     pressure: np.ndarray[shape[N, Z], np.dtype[float_]]
...     temperature: np.ndarray[shape[N, Z], np.dtype[float_]]
...     lcl_index: tuple[
...         np.ndarray[shape[N], np.dtype[np.intp]],
...         np.ndarray[shape[Z], np.dtype[np.intp]],
...     ]
...     @property
...     def lcl_pressure(self) -> np.ndarray[shape[N], np.dtype[float_]]:
...         return self.pressure[self.lcl_index]
...     @property
...     def lcl_temperature(self) -> np.ndarray[shape[N], np.dtype[float_]]:
...         return self.temperature[self.lcl_index]
...     def below_lcl(self, copy: bool = True):
...         P, T = self.pressure, self.temperature
...         N, Z = self.shape
...         if copy:
...             P, T = P.copy(), T.copy()
...         indices = np.ones(N, dtype=np.int_)[:, None] * np.arange(Z)
...         mask = indices > self.lcl_index[-1][:, None]
...         P[mask] = np.nan
...         T[mask] = np.nan
...         return P, T
...     def above_lcl(self, copy: bool = True):
...         P, T = self.pressure, self.temperature
...         N, Z = self.shape
...         if copy:
...             P, T = P.copy(), T.copy()
...         indices = np.ones(N, dtype=np.int_)[:, None] * np.arange(Z)
...         mask = indices < self.lcl_index[-1][:, None]
...         P[mask] = np.nan
...         T[mask] = np.nan
...         return P, T
...     @property
...     def shape(self) -> shape[N, Z]:
...         return self.pressure.shape  # type: ignore
... 
>>> def parcel_profile(
...     pressure: Annotated[Pascal[np.ndarray[shape[Z], np.dtype[float_]]], "pressure levels"],
...     temperature: Annotated[Kelvin[np.ndarray[shape[N], np.dtype[np.float_]]], "surface temperature"],
...     dewpoint: Annotated[Kelvin[np.ndarray[shape[N], np.dtype[np.float_]]], "surface dewpoint temperature"],
...     sfc_pressure: Annotated[Pascal[np.ndarray[shape[N], np.dtype[np.float_]]], "surface pressure"] | None = None,
... ) -> ParcelProfile[float_]:
...     # add a nan value to the end of the pressure array
...     dtype = pressure.dtype
...     pressure = np.append(pressure, np.nan)
...     N, Z = temperature.shape[0], pressure.shape[0]
...     p0 = sfc_pressure or pressure[:1].repeat(N)
...     assert temperature.shape == dewpoint.shape == p0.shape == (N,)
...     indices = np.arange(N)
...     # - calculate LCL
...     lcl_p, lcl_t = nzt.lcl(p0, temperature, dewpoint)  # (N,)
...     mask = pressure >= lcl_p.reshape(-1, 1)  # (N, Z)
...     mask[indices, np.argmin(mask, axis=1) + 1] = 0
...     # [ pressure ]
...     P = np.full((N, Z), np.nan, dtype=dtype)
...     nx, zx = np.nonzero(mask)
...     P[nx, zx] = pressure[zx]
...     nx, zx = np.nonzero(~mask & ~np.isnan(pressure))
...     P[nx, zx + 1] = pressure[zx]
...     # [ temperature ]
...     T = np.where(
...         mask,
...         nzt.dry_lapse(np.where(mask, P, np.nan), temperature[:, np.newaxis], p0[:, np.newaxis]),  # lower
...         nzt.moist_lapse(np.where(~mask, P, np.nan), temperature, lcl_p),  # upper
...     )
...     lcl_index = np.nonzero(np.isnan(P))
...     assert len(lcl_index) == 2 and len(lcl_index[0]) == N == len(lcl_index[1])
...     P[lcl_index] = lcl_p
...     T[lcl_index] = lcl_t
...     return ParcelProfile(P, T, lcl_index)
... 
>>> pressure = np.array(
...     [1013, 1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 725, 700, 650, 600, 550, 500, 450, 400, 350, 300],
... ).astype(float)
>>> 
>>> pressure *= 100.0
>>> temperature = np.array(
...     [
...         [243, 242, 241, 240, 239, 237, 236, 235, 233, 232, 231, 229, 228, 226, 235, 236, 234, 231, 226, 221, 217, 211],
...         [250, 249, 248, 247, 246, 244, 243, 242, 240, 239, 238, 236, 235, 233, 240, 239, 236, 232, 227, 223, 217, 211],
...         [293, 292, 290, 288, 287, 285, 284, 282, 281, 279, 279, 280, 279, 278, 275, 270, 268, 264, 260, 254, 246, 237],
...         [300, 299, 297, 295, 293, 291, 292, 291, 291, 289, 288, 286, 285, 285, 281, 278, 273, 268, 264, 258, 251, 242],
...     ]
... ).astype(float)
>>> 
>>> dewpoint = np.array(
...     [
...         [224, 224, 224, 224, 224, 223, 223, 223, 223, 222, 222, 222, 221, 221, 233, 233, 231, 228, 223, 218, 213, 207],
...         [233, 233, 232, 232, 232, 232, 231, 231, 231, 231, 230, 230, 230, 229, 237, 236, 233, 229, 223, 219, 213, 207],
...         [288, 288, 287, 286, 281, 280, 279, 277, 276, 275, 270, 258, 244, 247, 243, 254, 262, 248, 229, 232, 229, 224],
...         [294, 294, 293, 292, 291, 289, 285, 282, 280, 280, 281, 281, 278, 274, 273, 269, 259, 246, 240, 241, 226, 219],
...     ]
... ).astype(float)
>>> 
>>> pp = parcel_profile(pressure, temperature[:, 0], dewpoint[:, 0])
>>> pp
ParcelProfile(pressure=array([[101300.        , 100000.        ,  97500.        ,
         95000.        ,  92500.        ,  90000.        ,
         87500.        ,  85000.        ,  82500.        ,
         80000.        ,  77500.        ,  75000.        ,
         72840.34476965,  72500.        ,  70000.        ,
         65000.        ,  60000.        ,  55000.        ,
         50000.        ,  45000.        ,  40000.        ,
         35000.        ,  30000.        ],
       [101300.        , 100000.        ,  97500.        ,
         95000.        ,  92500.        ,  90000.        ,
         87500.        ,  85000.        ,  82500.        ,
         80000.        ,  77500.        ,  75973.87549461,
         75000.        ,  72500.        ,  70000.        ,
         65000.        ,  60000.        ,  55000.        ,
         50000.        ,  45000.        ,  40000.        ,
         35000.        ,  30000.        ],
       [101300.        , 100000.        ,  97500.        ,
         95000.        ,  94064.316261  ,  92500.        ,
         90000.        ,  87500.        ,  85000.        ,
         82500.        ,  80000.        ,  77500.        ,
         75000.        ,  72500.        ,  70000.        ,
         65000.        ,  60000.        ,  55000.        ,
         50000.        ,  45000.        ,  40000.        ,
         35000.        ,  30000.        ],
       [101300.        , 100000.        ,  97500.        ,
         95000.        ,  92815.92020834,  92500.        ,
         90000.        ,  87500.        ,  85000.        ,
         82500.        ,  80000.        ,  77500.        ,
         75000.        ,  72500.        ,  70000.        ,
         65000.        ,  60000.        ,  55000.        ,
         50000.        ,  45000.        ,  40000.        ,
         35000.        ,  30000.        ]]), temperature=array([[243.        , 242.10489757, 240.35991212, 238.58266792,
        236.7716956 , 234.92541688, 233.04213329, 231.12001323,
        229.15707727, 227.15118134, 225.09999737, 223.00099101,
        221.1471183 , 242.69930395, 240.44395172, 235.6779225 ,
        230.55000379, 225.02949201, 219.07631918, 212.63555153,
        205.62970772, 197.94737236, 189.42470981],
       [250.        , 249.07911272, 247.2838602 , 245.45541967,
        243.59227942, 241.69281572, 239.75528117, 237.77779139,
        235.75830995, 233.69463101, 231.58435943, 230.27213523,
        249.18901695, 247.04805761, 244.81931713, 240.08476438,
        234.95691952, 229.40404586, 223.38706435, 216.85347517,
        209.72873003, 201.90376347, 193.21544608],
       [293.        , 291.92072011, 289.81668416, 287.67375186,
        286.85549918, 292.40477687, 291.42720671, 290.41517912,
        289.36607813, 288.2769818 , 287.14461382, 285.96528574,
        284.73482722, 283.44850185, 282.10090464, 279.19621118,
        275.9585755 , 272.30349147, 268.11347696, 263.22214704,
        257.39203344, 250.29057047, 241.48515066],
       [300.        , 298.89493527, 296.74063224, 294.54650361,
        292.58698851, 299.89311651, 299.03285565, 298.14571362,
        297.22985171, 296.28323374, 295.30359666, 294.2884156 ,
        293.23486191, 292.13975254, 290.99948846, 288.56658581,
        285.89567308, 282.93245038, 279.60178469, 275.79587358,
        271.3531953 , 266.01916996, 259.36993871]]), lcl_index=(array([0, 1, 2, 3]), array([12, 11,  4,  4])))
>>> pp.lcl_temperature
array([221.1471183 , 230.27213523, 286.85549918, 292.58698851])
>>> pp.lcl_pressure
array([72840.34476965, 75973.87549461, 94064.316261  , 92815.92020834])
>>> pp.above_lcl()
(array([[           nan,            nan,            nan,            nan,
                   nan,            nan,            nan,            nan,
                   nan,            nan,            nan,            nan,
        72840.34476965, 72500.        , 70000.        , 65000.        ,
        60000.        , 55000.        , 50000.        , 45000.        ,
        40000.        , 35000.        , 30000.        ],
       [           nan,            nan,            nan,            nan,
                   nan,            nan,            nan,            nan,
                   nan,            nan,            nan, 75973.87549461,
        75000.        , 72500.        , 70000.        , 65000.        ,
        60000.        , 55000.        , 50000.        , 45000.        ,
        40000.        , 35000.        , 30000.        ],
       [           nan,            nan,            nan,            nan,
        94064.316261  , 92500.        , 90000.        , 87500.        ,
        85000.        , 82500.        , 80000.        , 77500.        ,
        75000.        , 72500.        , 70000.        , 65000.        ,
        60000.        , 55000.        , 50000.        , 45000.        ,
        40000.        , 35000.        , 30000.        ],
       [           nan,            nan,            nan,            nan,
        92815.92020834, 92500.        , 90000.        , 87500.        ,
        85000.        , 82500.        , 80000.        , 77500.        ,
        75000.        , 72500.        , 70000.        , 65000.        ,
        60000.        , 55000.        , 50000.        , 45000.        ,
        40000.        , 35000.        , 30000.        ]]), array([[         nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
        221.1471183 , 242.69930395, 240.44395172, 235.6779225 ,
        230.55000379, 225.02949201, 219.07631918, 212.63555153,
        205.62970772, 197.94737236, 189.42470981],
       [         nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan, 230.27213523,
        249.18901695, 247.04805761, 244.81931713, 240.08476438,
        234.95691952, 229.40404586, 223.38706435, 216.85347517,
        209.72873003, 201.90376347, 193.21544608],
       [         nan,          nan,          nan,          nan,
        286.85549918, 292.40477687, 291.42720671, 290.41517912,
        289.36607813, 288.2769818 , 287.14461382, 285.96528574,
        284.73482722, 283.44850185, 282.10090464, 279.19621118,
        275.9585755 , 272.30349147, 268.11347696, 263.22214704,
        257.39203344, 250.29057047, 241.48515066],
       [         nan,          nan,          nan,          nan,
        292.58698851, 299.89311651, 299.03285565, 298.14571362,
        297.22985171, 296.28323374, 295.30359666, 294.2884156 ,
        293.23486191, 292.13975254, 290.99948846, 288.56658581,
        285.89567308, 282.93245038, 279.60178469, 275.79587358,
        271.3531953 , 266.01916996, 259.36993871]]))
>>> pp.below_lcl()
(array([[101300.        , 100000.        ,  97500.        ,
         95000.        ,  92500.        ,  90000.        ,
         87500.        ,  85000.        ,  82500.        ,
         80000.        ,  77500.        ,  75000.        ,
         72840.34476965,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan],
       [101300.        , 100000.        ,  97500.        ,
         95000.        ,  92500.        ,  90000.        ,
         87500.        ,  85000.        ,  82500.        ,
         80000.        ,  77500.        ,  75973.87549461,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan],
       [101300.        , 100000.        ,  97500.        ,
         95000.        ,  94064.316261  ,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan],
       [101300.        , 100000.        ,  97500.        ,
         95000.        ,  92815.92020834,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan,             nan,
                    nan,             nan]]), array([[243.        , 242.10489757, 240.35991212, 238.58266792,
        236.7716956 , 234.92541688, 233.04213329, 231.12001323,
        229.15707727, 227.15118134, 225.09999737, 223.00099101,
        221.1471183 ,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan],
       [250.        , 249.07911272, 247.2838602 , 245.45541967,
        243.59227942, 241.69281572, 239.75528117, 237.77779139,
        235.75830995, 233.69463101, 231.58435943, 230.27213523,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan],
       [293.        , 291.92072011, 289.81668416, 287.67375186,
        286.85549918,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan],
       [300.        , 298.89493527, 296.74063224, 294.54650361,
        292.58698851,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan,          nan]]))
>>>