search

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.25k stars 416 forks source link

metpy.calc.thermo mask and data must be the same size #972

Closed am-thyst closed 5 years ago

am-thyst commented 5 years ago 
mb900 = xr.open_dataset("mb900_data.nc")
sfc_1 = xr.open_dataset("sfc_data.nc")
sfc_2 = xr.open_dataset("sfc_data2.nc")
dp_900 = np.load("/nerc/n02/n02/amethyst/npys/Td.npy") #900 hPa dewpoint (calculated using atmos)

####setting the variables
sfc_p = np.array(sfc_1.sp) #surface pressure
sfc_p = sfc_p / 100 #Pa to hPa
#create an array of same dimensions all set to 900hPa
pressure_900 = sfc_p.copy()
pressure_900[:,:] = 900

sfc_t = np.array(sfc_2.t2m) #temperature at surface
sfc_dp = np.array(sfc_2.d2m) #dewpoint at surface 
temp_900 = np.array(mb900.t) #teperature at 900hPa

####combining sfc pressure values and 900 hPa values for the parcel_profile function
pp_pressure = np.array([sfc_p, pressure_900]) * units.hPa
pp_temp = np.array([sfc_t, temp_900]) * units.kelvin
pp_dp = np.array([sfc_dp, dp_900]) * units.kelvin

####
parcel_profile = parcel_profile(pp_pressure, pp_temp, pp_dp)

/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/metpy/calc/thermo.py:690: RuntimeWarning: invalid value encountered in log
  val = np.log(e / sat_pressure_0c)
/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/pint/quantity.py:1403: RuntimeWarning: invalid value encountered in log
  out = uf(*mobjs)
Traceback (most recent call last):
  File "/var/spool/torque/mom_priv/jobs/65276.rdf-xcat.SC", line 45, in <module>
    parcel_profile = parcel_profile(pp_pressure, pp_temp, pp_dp)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/metpy/xarray.py", line 381, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/metpy/calc/thermo.py", line 466, in parcel_profile
    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpt)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/metpy/calc/thermo.py", line 527, in _parcel_profile_helper
    press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpt)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/metpy/xarray.py", line 381, in wrapper
    return func(*args, **kwargs)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/metpy/calc/thermo.py", line 315, in lcl
    xtol=eps, maxiter=max_iters)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 904, in fixed_point
    return _fixed_point_helper(func, x0, args, xtol, maxiter, use_accel)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 853, in _fixed_point_helper
    relerr = _lazywhere(p0 != 0, (p, p0), f=_relerr, fillvalue=p)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/scipy/_lib/_util.py", line 54, in _lazywhere
    np.place(out, cond, f(*temp))
  File "/nerc/n02/n02/amethyst/miniconda3/envs/myenv/lib/python2.7/site-packages/numpy/lib/function_base.py", line 1607, in place
    return _insert(arr, mask, vals)
ValueError: place: mask and data must be the same size

Shape for all of the components is the same (28490, 241, 480) Could someone let me know what I'm doing wrong?

am-thyst commented 5 years ago

ping @dopplershift

jrleeman commented 5 years ago

@am-thyst - sorry we're a bit swamped getting ready for AMS. I'll try to dig into this later today.

jrleeman commented 5 years ago

@am-thyst so we currently don't work on grids - you'll have to manually loop over the grid points which we admit will be rather slow. We'd like to speed up grid calculations, and hopefully XArray can help us moving forward.

am-thyst commented 5 years ago

@jrleeman ah, that's a sticky one. I've been converting them all to numpy arrays so far, I think I read somewhere that metpy doesn't work with xarray? But for sure in the future this would be a good change.

I think the cape calculation works like this in most python packages. Forget parcel_profile() for a second, let's take a different example of surface_based_cape_cin(), I got a similar error for probably the same reason. I definitely have the computing power to go over every grid point, so please would you be able to give me a few pointers on how to approach this (with sbcape/cin), inputting every grid point back into an array of the original shape (28490, 241, 480)?

jrleeman commented 5 years ago

We do use XArray, just not everywhere yet. The transition is happening, but there are some pesky unit issues and other things that need to settle out in the ecosystem before it's ready to go as a full on MetPy data model. You might have a look at our training materials, which include an xarray module and it is used a few places throughout (satellite being one example). https://github.com/Unidata/python-workshop

Here's a looping example (untested) that assumes your data are x, y, p:

rows, columns, _ = np.shape(my_data)
results = np.zeros((rows, columns))

for i in range(rows):
    for j in range(columns):
         profile = my_data[i, j, :]
         # Do calculations, whatever
         profile_result = ......
         results[i, j] = profile_result
am-thyst commented 5 years ago

@jrleeman thanks for your response, I thought I better test this out so I understand it and to make sure it does what I want it to do before applying it to my huge dataset. I just made up some values and tried to make it somewhat similar to what I'm doing already, just more simple.

Would you mind having a quick look? I think it's the last line I have a problem with, giving ValueError: setting an array element with a sequence, it's just I'm unsure how to put the values into their respective positions in the array. If you see I've done anything really wrong then please let me know!

sfc_p = np.array([1013, 1010, 1014])
p_900 = np.array([900, 900, 900])
sfc_t = np.array([285, 286, 284])
t_900 = np.array([275, 273, 278])
sfc_dp = np.array([280, 279, 281])
dp_900 = np.array([270, 268, 273])

p = np.array([sfc_p, p_900]) * units.hPa
t = np.array([sfc_t, t_900]) * units.kelvin
dp = np.array([sfc_dp, dp_900]) * units.kelvin

_, value = np.shape(t)
results = np.zeros(value)

for i in range(value):
    temp = t[:,i]
    pressure = p[:,i]
    dewpoint = dp[:,i]
    cape, cin = surface_based_cape_cin(pressure, temp, dewpoint)
    results[i] = cape
jrleeman commented 5 years ago

In this example your arrays are all one dimensional and you're indexing into them as if they were two dimensional. Have a look at our NumPy training as well as the numpy quickstart

To give you an idea, let's say your data are in x, y, pressure order... temperatures might look like:

temperatures = np.array([[[22, 23, 24], [24, 25, 26]],
[[27, 28, 29], [30, 32, 34]]])
am-thyst commented 5 years ago

@jrleeman I'm indexing them as two dimensional because I've combined the surface values with the pressure values, eg:

p = np.array([sfc_p, p_900]) * units.hPa #1D sfc_p array combined with 1D 900 mb pressure array

print(p)
<Quantity([[1013 1010 1014]
 [ 900  900  900]], 'hectopascal')>

making their shape 2 dimensional (2, 3)

jrleeman commented 5 years ago

Either way - you're going to use the indexing strategy I've shown to iterate across your data grids.

am-thyst commented 5 years ago

@jrleeman I don't really get it, what I interpret from your code example is that the data will have rows, columns (essentially longitudes and latitudes?) and another dimension which isn't represented in the final results? Whereas my data has another dimension - time, so here's what I guessed your code would translate to when applied to my data:

_, time, latitude, longitude = np.shape(temperature) # shape (2, 28490, 241, 480)
results = np.zeros((time, latitude, longitude))

for i in range(time):
    for j in range(latitude):
        for k in range(longitude):
            temperature[:, i, j, k] # same for the other two inputs
            cape, cin = surface_based_cape_cin(pressure, temperature, dewpoint)
            results[i, j, k] = cape

Am I missing something? Sorry for dragging this out

jrleeman commented 5 years ago

surface_based_cape_cin takes pressure, temperature, dewpoint profiles, so you need a height/pressure dimension.

https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.surface_based_cape_cin.html#metpy.calc.surface_based_cape_cin

am-thyst commented 5 years ago

@jrleeman the pressure dimension is the 2 in (2, 28490, 241, 480), so temp[0] would be values for the surface and temp[1] would be values for the temperature at 900hpa

jrleeman commented 5 years ago

So you only have two pressures? I doubt you'll be able to be meaningful results from that. You can use np.squeeze if you need to remove zero size dimensions.

am-thyst commented 5 years ago

@jrleeman I have all pressure levels available to me, so I can definitely add more to the pressure dimension to improve the results, it's just the loop code I'm kind of confused with at the minute

jrleeman commented 5 years ago

Assign the slice you took to a variable:

_, time, latitude, longitude = np.shape(temperature) # shape (2, 28490, 241, 480)
results = np.zeros((time, latitude, longitude))

for i in range(time):
    for j in range(latitude):
        for k in range(longitude):
            temperature_profile = temperature[:, i, j, k] # same for the other two inputs
            cape, cin = surface_based_cape_cin(pressure, temperature_profile, dewpoint_profile)
            results[i, j, k] = cape
am-thyst commented 5 years ago

@jrleeman thank you I was almost there! I'll try it tomorrow 👍

am-thyst commented 5 years ago

@jrleeman 

/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/interpolate/one_dimension.py:148: UserWarning: Interpolation point out of data bounds encountered
  warnings.warn('Interpolation point out of data bounds encountered')
/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/pint/quantity.py:1070: RuntimeWarning: invalid value encountered in less
  return op(self._magnitude, other._magnitude)
/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/calc/thermo.py:1406: RuntimeWarning: invalid value encountered in greater
  keep_idx = np.ediff1d(x, to_end=[1]) > 0
/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/pint/quantity.py:1065: RuntimeWarning: invalid value encountered in less
  return op(self._convert_magnitude_not_inplace(UnitsContainer()), other)
/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/pint/quantity.py:1065: RuntimeWarning: invalid value encountered in greater
  return op(self._convert_magnitude_not_inplace(UnitsContainer()), other)
/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/interpolate/one_dimension.py:138: UserWarning: Interpolation point out of data bounds encountered
  warnings.warn('Interpolation point out of data bounds encountered')
Traceback (most recent call last):
  File "/var/spool/torque/mom_priv/jobs/65714.rdf-xcat.SC", line 54, in <module>
    cape, cin = surface_based_cape_cin(pressure_profile, temperature_profile, dewpoint_profile)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/xarray.py", line 381, 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 1653, 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 381, 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 1332, 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 381, 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 374, in lfc
    x, y = lcl(pressure[0], temperature[0], dewpt[0])
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/metpy/xarray.py", line 381, 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 315, in lcl
    xtol=eps, maxiter=max_iters)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 904, in fixed_point
    return _fixed_point_helper(func, x0, args, xtol, maxiter, use_accel)
  File "/nerc/n02/n02/amethyst/miniconda3/envs/py3env/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 858, in _fixed_point_helper
    raise RuntimeError(msg)
RuntimeError: Failed to converge after 50 iterations, value is nan
jrleeman commented 5 years ago

This is indicating that a calculation isn't converging - likely because of something in your data. I'd advise you to pull a couple of profiles from the dataset (maybe even the problematic one) and plot them on a skewT to make sure that you are getting what you think you are.

am-thyst commented 5 years ago

I believe this was solved by adding data from more pressure levels. Will update if problem reoccurs.