Stream: python-questions

Topic: usage interp_hybrid_to_pressure

view this post on Zulip Brian Medeiros (May 04 2021 at 19:29):

Hi all, but especially @Orhan Eroglu :
I'm getting an error numpy.AxisError: axis 1 is out of bounds for array of dimension 1 when trying to use GeoCAT's interp_hybrid_to_pressure. Here's what I'm trying:

import geocat.comp as gc
import xarray as xr
loc = "/glade/p/cesm/pcwg/jenkay/COSP/cesm21/f.e21.FHIST_BGC.f09_f09_mg17.CMIP6-AMIP.001_cosp1/atm/proc/tseries/month_1/"
ds_ps = xr.open_dataset(loc+"f.e21.FHIST_BGC.f09_f09_mg17.CMIP6-AMIP.001_cosp1.cam.h0.PS.197901-201412.nc")
ds_t = xr.open_dataset(loc+"f.e21.FHIST_BGC.f09_f09_mg17.CMIP6-AMIP.001_cosp1.cam.h0.T.197901-201412.nc")
temp = ds_t["T"]
ps = ds_ps["PS"]
temp_plev = gc.interp_hybrid_to_pressure(temp, ps, ds_t['hyam'], ds_t['hybm'], lev_dim='lev')

This is with python 3.9.2 in a brand new geocat environment on casper (with the addition of conda install netcdf4).

Any idea what is happening?

view this post on Zulip Deepak Cherian (May 04 2021 at 19:33):

https://github.com/NCAR/geocat-comp/blob/eae0d27f353db2bfae9ae36ad2b615fdca66cf44/src/geocat/comp/interp_hybrid_to_pressure.py#L70

This is a bug. interp_axis should always be -1, since lev_dim is a core dim, apply_ufunc will transpose so that lev_dim is the last dimension.

@Brian Medeiros try gc.interp_hybrid_to_pressure(temp.transpose(..., "lev"), ...).

view this post on Zulip Brian Medeiros (May 04 2021 at 21:38):

Thanks, @Deepak Cherian ! That makes sense. Applying transpose appears to be working.

view this post on Zulip Brian Medeiros (May 04 2021 at 22:22):

I spoke too soon. I can get it to work for a single column, but it seems like as soon as I use an actual multidimensional array it is giving a similar kind of error. I tried to use a subset of the data:

ntime = 2
temp_sample = temp.isel(time=slice(None,ntime)).transpose("time","lat","lon","lev")
ps_sample = ps.isel(time=slice(None,ntime))
print(temp_sample.shape)  # (2, 192, 288, 32)
print(ps_sample.shape)  # (2, 192, 288)
tmp = gc.interp_hybrid_to_pressure(temp_sample, ps_sample, ds_t['hyam'], ds_t['hybm'])
tmp

and it says AxisError: axis 3 is out of bounds for array of dimension 1

view this post on Zulip Brian Medeiros (May 05 2021 at 02:52):

Well, I don't completely understand what is going on, but I can share a workaround for this issue. I switched from the MetPy interpolation to just plain old numpy interp:

    if method == 'linear':
#         func_interpolate = metpy.interpolate.interpolate_1d
        func_interpolate = np.interp

and corresponding adjustment of the vertical remapper:

    def _vertical_remap(data, pressure):
        """Define interpolation function."""
        return func_interpolate(new_levels, pressure, data, left=np.nan, right=np.nan)

This appears to behave nicely.

view this post on Zulip Orhan Eroglu (May 05 2021 at 15:36):

I will be looking into it soon

view this post on Zulip Ryan May (May 11 2021 at 19:53):

If you find we're doing something wrong, please don't hesitate to open an issue.

view this post on Zulip Orhan Eroglu (May 18 2021 at 21:35):

Especially @Brian Medeiros and @Ryan May , and also others, sorry for the late report back. I have had some findings now:

1) The pressure and data were not of the same shape (after _pressure_from_hybrid call) though it is required by metpy.interpolate.interpolate_1d. I first corrected this one:

# Calculate pressure levels at the hybrid levels
    pressure = _pressure_from_hybrid(ps, hyam, hybm, p0)  # Pa

    pressure = pressure.transpose(*data.dims)

2) After that, the code was still failing within the xarray.apply_ufunc call. I wanted to first see how metpy.interpolate.interpolate_1d or metpy.interpolate.log_interpolate_1d is behaving when data is of type numpy.ndarray or xarray.DataArray. So, I tested the following cases:

a- Numpy Inputs: When I converted Brian's inputs to Numpy and called metpy functions directly (i.e. no xarray.apply_ufunc), the code ran successfully.

    # Calling with Numpy is successful
    output_np = _vertical_remap(data.values, pressure.values)

b- Xarray Inputs: When I kept Brian's inputs as is ( i.e. Xarray) and called metpy functions directly (i.e. no xarray.apply_ufunc), the code failed.

    # Calling with Xarray is failing!!!
    output_xr = _vertical_remap(data, pressure)

Here is the failure message, which is clearer than the failure messages that were output while calling via xarray.apply_ufunc (Though, the original failure Brian found in the original code, i.e. via xarray.apply_ufunc, was due to this same xarray thing I think):

  File "/Users/oero/miniconda3/envs/geocat_comp_build/lib/python3.9/site-packages/geocat/comp/interp_hybrid_to_pressure.py", line 89, in _vertical_remap
    return func_interpolate(new_levels, pressure, data, axis=interp_axis)
  File "/Users/oero/miniconda3/envs/geocat_comp_build/lib/python3.9/site-packages/metpy/xarray.py", line 1206, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "/Users/oero/miniconda3/envs/geocat_comp_build/lib/python3.9/site-packages/metpy/interpolate/one_dimension.py", line 127, in interpolate_1d
    minv = np.apply_along_axis(np.searchsorted, axis, xp, x[sort_x])
  File "<__array_function__ internals>", line 5, in apply_along_axis
TypeError: no implementation found for 'numpy.apply_along_axis' on types that implement __array_function__: [<class 'pint.quantity.build_quantity_class.<locals>.Quantity'>]

I couldn't dig into the actual issue with metpy & xarray for the time; I will continue my investigation into it.

@Ryan May : I can create an issue in metpy repo, but I am not sure how to provide them/you with an easily replicable GeoCAT-comp run of this (i.e. it may require you to setup GeoCAT-comp etc.). Maybe I could create a branch/PR from these changes I have made under our code, and metpy folks would like to step in? Let me know what is more convenient for you (if you think it is a metpy issue).

view this post on Zulip Orhan Eroglu (May 27 2021 at 03:21):

Hi @Brian Medeiros

Please see GeoCAT-comp PR #155 for a workaround to get interp_hybrid_to_pressure successfully for Xarray inputs of 4D or more dims. With this new version of the code, your original script (i.e. without transposing data) should work successfully. I already tested it. Please let me know if you can test it with the code in the PR; otherwise, you can also give it a try with our newer version via conda around end of May or early June. By the way, I tested both unchunked and chunked versions of your data and got dramatically different performances, which I could share with you.

@Ryan May I think that this is related to metpy.interpolate's xarray interface; thus, I have created an issue with a minimal code that does directly run on metpy functions (i.e. no geocat-comp functionality). Please see the Metpy Issue #1889 here

Last updated: Jan 30 2022 at 12:01 UTC