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.24k stars 413 forks source link

index out of bounds sbcape #1011

Open am-thyst opened 5 years ago

am-thyst commented 5 years ago

Kind of follows on from issues #972 and #997 -

With a simple example, the surface_based_cape_cin function works fine:

from metpy.calc.thermo import surface_based_cape_cin
from metpy.units import units

t = [288, 284, 278, 272, 265, 260, 254, 250] * units.kelvin
dp = [285, 278, 269, 266, 261, 254, 250, 244] * units.kelvin
p = [1013, 950, 900, 800, 750, 700, 650, 600] * units.hPa

 surface_based_cape_cin(p, t, dp)
(<Quantity(925.6897878888357, 'joule / kilogram')>, <Quantity(-8.118571674744551, 'joule / kilogram')>)

but when applying it to my loop, where pressure, temperature & dewpoint are definitely of shape (8,):

dict = {}
# Calculates the timesteps between date 1 and date 2
datetime_range = datetime_calculation((datetime(1979, 1, 1)), (datetime(1979, 12, 31)))

arra = np.zeros((int(datetime_range), 241, 480))
# Produces an array of intended output shape which can be used for multiprocessing:
counter = sharedmem.empty_like(arra) 

sfc_1 = xr.open_dataset("sfc_data.nc")
sfc_2 = xr.open_dataset("sfc_data2.nc")
sfc_p = np.array(sfc_1.msl) / 100
sfc_t = np.array(sfc_2.t2m)
sfc_dp = np.array(sfc_2.d2m)

for x in range(650, 1000, 50): # Opens the xarray data arrays for each pressure level & parameter
  dict["p" + str(x)] = pressure_levels(int(x))
  dict["t" + str(x)] = temps(str(x))
  dict["dp" + str(x)] = dewpoints(str(x))

with sharedmem.MapReduce(np=59) as pool:
  def cape_func(abc):
    for j in range(int(abc)-1, int(abc)):
      for k in range(0, 241):  
        for l in range(0, 480):
          pressure = [sfc_p[j,k,l], dict["p950"][j,k,l], dict["p900"][j,k,l], dict["p850"][j,k,l], dict["p800"][j,k,l], dict["p750"][j,k,l], dict["p700"][j,k,l], dict["p650"][j,k,l]] * units.hPa
          temperature = [sfc_t[j,k,l], dict["t950"][j,k,l], dict["t900"][j,k,l], dict["t850"][j,k,l], dict["t800"][j,k,l], dict["t750"][j,k,l], dict["t700"][j,k,l], dict["t650"][j,k,l]] * units.kelvin
          dewpoint = [sfc_dp[j,k,l], dict["dp950"][j,k,l], dict["dp900"][j,k,l], dict["dp850"][j,k,l], dict["dp800"][j,k,l], dict["dp750"][j,k,l], dict["dp700"][j,k,l], dict["dp650"][j,k,l]] * units.kelvin
          cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint)
          cape = cape / (units.joule / units.kilogram)
          counter[j,k,l] = float(cape)
  pool.map(cape_func, range(1, (int(datetime_range)+1))

I get the following traceback:

Traceback (most recent call last):
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/sharedmem/sharedmem.py", line 294, in _slaveMain
    self.main(self, *self.args)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/sharedmem/sharedmem.py", line 628, in _main
    r = realfunc(work)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/sharedmem/sharedmem.py", line 703, in realfunc
    else: return func(i)
  File "/var/spool/torque/mom_priv/jobs/237175.rdf-xcat.SC", line 53, in cape_func
    cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/xarray.py", line 436, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/units.py", line 305, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/calc/thermo.py", line 1697, in surface_based_cape_cin
    return cape_cin(p, t, td, profile)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/xarray.py", line 436, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/units.py", line 305, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/calc/thermo.py", line 1374, in cape_cin
    parcel_temperature_profile=parcel_profile)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/xarray.py", line 436, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/units.py", line 305, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/calc/thermo.py", line 424, in lfc
    return x[0], y[0]
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/pint/quantity.py", line 1281, in __getitem__
    value = self._magnitude[key]
IndexError: index 0 is out of bounds for axis 0 with size 0

edit: works when taking just one point from the loop, so definitely not the data that's the problem?

j, k, l = 0, 200, 400
pressure = [sfc_p[j,k,l], dict["p950"][j,k,l], dict["p900"][j,k,l], dict["p850"][j,k,l], dict["p800"][j,k,l], dict["p750"][j,k,l], dict["p700"][j,k,l], dict["p650"][j,k,l]] * units.hPa
temperature = [sfc_t[j,k,l], dict["t950"][j,k,l], dict["t900"][j,k,l], dict["t850"][j,k,l], dict["t800"][j,k,l], dict["t750"][j,k,l], dict["t700"][j,k,l], dict["t650"][j,k,l]] * units.kelvin
dewpoint = [sfc_dp[j,k,l], dict["dp950"][j,k,l], dict["dp900"][j,k,l], dict["dp850"][j,k,l], dict["dp800"][j,k,l], dict["dp750"][j,k,l], dict["dp700"][j,k,l], dict["dp650"][j,k,l]] * units.kelvin
cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint)
print(cape)

11.299586587859002 joule / kilogram
am-thyst commented 5 years ago

ping @dopplershift

dopplershift commented 5 years ago

My guess is that one particular column from the grid has some conditions that our CAPE code is not happy with. Any way to identify which set of indices is being worked on when the error occurs?

sgdecker commented 5 years ago

I think this is the same as #945, for which I have a PR.

am-thyst commented 5 years ago

@sgdecker any way to work around it?

sgdecker commented 5 years ago

The fix will be in master soon. Since the problem is internal to the lfc function, the only thing I can suggest other than waiting is to put your call to surface_based_cape_cin in a try/except block so that at least the gridpoints that don't trigger the bug can be computed.

zbruick commented 5 years ago

@am-thyst did the PR provided by @sgdecker fix this issue? Just going back through issues and trying to figure out what hasn't been solved yet.

zbruick commented 5 years ago

Hearing nothing on this, I'm going to close this issue.

am-thyst commented 5 years ago

Hi @zbruick sorry for the late response. I'm still getting this error, even with the PR implemented. It also occurs on the most_unstable_cape_cin function.

zbruick commented 5 years ago

To confirm, you've updated to the current version of master and tested again? Just want to make sure before investigating more, as the PR hasn't been released yet, but will be in 0.11 when that's released shortly.

michaelsessa commented 5 years ago

I also get this error with all CAPE related calculations ie. most_unstable_cape_cin, surface_based_cape_cin.

zbruick commented 5 years ago

Hey @michaelsessa! Can you provide any more details and do you have a reproducible example for this? Is it erring out for any point, or are you doing this on a grid? If we can get a reproducible example for this problem, that would be a great help in solving it!

michaelsessa commented 5 years ago

Yes! I am calculating different environmental parameters across a domain of 72 grid points. It varies case to case with some giving me no errors at all and others successfully completing calculations over a portion of the grid points before producing the error. image

sgdecker commented 5 years ago

Given that line 424 of thermo.py is not a return statement in master, I will second https://github.com/Unidata/MetPy/issues/1011#issuecomment-532379979 I suspect this is 0.10, which doesn't include the PR that fixes this issue.

zbruick commented 5 years ago

@michaelsessa, thanks for the traceback. Could you test against the master branch of MetPy (install it in a new Conda environment)? If you need instructions, you can find them here. There was a PR that might have fixed this issue, and any testing that could be done to verify that would be very helpful! Otherwise, feel free to drop me your data file via email and I can play with it and see what's up.

michaelsessa commented 5 years ago

I created a new Conda environment and tested against the master branch of Metpy and the issue appears to be resolved!

zbruick commented 5 years ago

Ok thanks for checking...I'm glad to hear it!

dopplershift commented 4 years ago

@am-thyst Given that #1172 is merged (fixing #1167), is your original issue here still present on master?

dopplershift commented 4 years ago

Bumping from 0.11, but please do let me know if this still a problem.

nmjurg commented 4 years ago

Sorry for re-upping an old thread, but I get an "index out of bounds" error when attempting to calculate wet bulb temperatures. Oddly, it only occurs with certain sounding data (i.e. certain dates/observations):

works for:

doesn't work for:

Dates are generally selected at random except for a few local severe events. Any ideas?

============================================================ Code snippet:

df = WyomingUpperAir.request_data(date, station) #pressure, height, temp, dp, direction, speed, u_wind, v_wind, station, station_number, time, lat, lon, elev

print(df)

p =             df['pressure'].values * units.hPa
T =             df['temperature'].values * units.degC
Td =            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)
mask =          p >= 100 * units.hPa

Tw = mpcalc.wet_bulb_temperature(p, T, Td).to('degC')

============================================================ Error:

Traceback (most recent call last):
  File ".\MetPy_Sounding.py", line 42, in <module>
    Tw = mpcalc.wet_bulb_temperature(p, T, Td).to('degC')
  File "C:\ProgramData\Anaconda3\lib\site-packages\metpy\xarray.py", line 655, in wrapper
    return func(*args, **kwargs)
  File "C:\ProgramData\Anaconda3\lib\site-packages\metpy\units.py", line 320, in wrapper
    return func(*args, **kwargs)
  File "C:\ProgramData\Anaconda3\lib\site-packages\metpy\calc\thermo.py", line 2413, in wet_bulb_temperature
    ret[...] = moist_adiabat_temperatures[-1].magnitude
  File "C:\ProgramData\Anaconda3\lib\site-packages\pint\quantity.py", line 1719, in __getitem__
    return type(self)(self._magnitude[key], self._units)
IndexError: index -1 is out of bounds for axis 0 with size 0
sgdecker commented 4 years ago

Shouldn't this be a new issue? The original problem here was with sbcape not wet_bulb_temperature.

dopplershift commented 4 years ago

I went ahead and copied that comment to #1332.