Reference pressure array in geostrophy.geo_strf_dyn_height

[ledm]

The documentation for geo_strf_dyn_height says that p_ref can be a "float or array-like". However, if this function is given a numpy array, it leads to the error:

""" File "/users/modellers/ledm/anaconda3/envs/ar6/lib/python3.8/site-packages/gsw/geostrophy.py", line 66, in geo_strf_dyn_height p_ref = np.float(p_ref) TypeError: only size-1 arrays can be converted to Python scalars """

This is due to the cast to float on L66, which can not be applied to an array, a list, or a tuple: https://github.com/TEOS-10/GSW-Python/blob/master/gsw/geostrophy.py#L66

The simplest solution would be to change the documentation, but having an array-like input for p_ref would be nice.

[ledm]

This is my solution. It basically checks whether p_ref is an array, and if so, creates a float with p_ref.

Hope this helps!

""" def geo_strf_dyn_height(SA, CT, p, p_ref=0, axis=0, max_dp=1.0, interp_method='pchip'): """ Dynamic height anomaly as a function of pressure.

Parameters
----------
SA : array-like
    Absolute Salinity, g/kg
CT : array-like
    Conservative Temperature (ITS-90), degrees C
p : array-like
    Sea pressure (absolute pressure minus 10.1325 dbar), dbar
p_ref : float or array-like, optional
    Reference pressure, dbar
axis : int, optional, default is 0
    The index of the pressure dimension in SA and CT.
max_dp : float
    If any pressure interval in the input p exceeds max_dp, the dynamic
    height will be calculated after interpolating to a grid with this
    spacing.
interp_method : string {'pchip', 'linear'}
    Interpolation algorithm.

Returns
-------
dynamic_height : array
    This is the integral of specific volume anomaly with respect
    to pressure, from each pressure in p to the specified
    reference pressure.  It is the geostrophic streamfunction
    in an isobaric surface, relative to the reference surface.

"""
interp_methods = {'pchip' : 2, 'linear' : 1}
if interp_method not in interp_methods:
    raise ValueError('interp_method must be one of %s'
                     % (interp_methods.keys(),))
if SA.shape != CT.shape:
    raise ValueError('Shapes of SA and CT must match; found %s and %s'
                     % (SA.shape, CT.shape))
if p.ndim == 1 and SA.ndim > 1:
    if len(p) != SA.shape[axis]:
        raise ValueError('With 1-D p, len(p) must be SA.shape[axis];\n'
                         ' found %d versus %d on specified axis, %d'
                         % (len(p), SA.shape[axis], axis))
    ind = [np.newaxis] * SA.ndim
    ind[axis] = slice(None)
    p = p[tuple(ind)]
if isinstance(p_ref, int) or isinstance(p_ref, str) :
    p_ref = np.float(p_ref)
with np.errstate(invalid='ignore'):
    # The need for this context seems to be a bug in np.ma.any.
    if np.ma.any(np.ma.diff(np.ma.masked_invalid(p), axis=axis) <= 0):
        raise ValueError('p must be increasing along the specified axis')
p = np.broadcast_to(p, SA.shape)
goodmask = ~(np.isnan(SA) | np.isnan(CT) | np.isnan(p))
dh = np.empty(SA.shape, dtype=float)
dh.fill(np.nan)

order = 'F' if SA.flags.fortran else 'C'
for ind in indexer(SA.shape, axis, order=order):
    igood = goodmask[ind]
    # If p_ref is below the deepest value, skip the profile.
    pgood = p[ind][igood]
    if isinstance(p_ref, float):
        prefv = p_ref
    else:
        prefv = p_ref[ind[1]]
    if  len(pgood) > 1 and pgood[-1] >= prefv:
        sa = SA[ind][igood]
        ct = CT[ind][igood]
        # Temporarily add a top (typically surface) point and mixed layer
        # if p_ref is above the shallowest pressure.
        if pgood[0] > prefv:
            ptop = np.arange(prefv, pgood[0], max_dp)
            ntop = len(ptop)
            sa = np.hstack(([sa[0]]*ntop, sa))
            ct = np.hstack(([ct[0]]*ntop, ct))
            pgood = np.hstack((ptop, pgood))
        else:
            ntop = 0
        dh_all = _gsw_ufuncs.geo_strf_dyn_height_1(
                                     sa, ct, pgood, prefv, max_dp,
                                     interp_methods[interp_method])
        if ntop > 0:
            dh[ind][igood] = dh_all[ntop:]
        else:
            dh[ind][igood] = dh_all

return dh

"""

[efiring]

Is there an important use case for p_ref to vary from profile to profile? If not, I'm inclined to just fix the documentation. Supporting p_ref as an array is possible but it is likely to be confusing because of the generality of the inputs, in which an array p_ref would have one dimension fewer than SA etc., and the omitted dimension is given by 'axis'. I'm not sure the complexity would be worthwhile. Also, the calculation is being done with an internal python loop over the profiles, so if the user wants different p_ref values for different profiles, there is little penalty to doing the loop externally.

[efiring]

prefv = p_ref[ind[1]]

I think you are making an assumption here about the dimensions of the input and the value of 'axis'.

[ledm]

The use case for us is that if we have an array of profiles with varying sea floor depth, and we want the reference depth to be either the seafloor or 2000m, whichever is shallower.

[efiring]

I think this is a misunderstanding of the way the function works, and how to use it. p_ref is the lowest-pressure limit of the integral, not the highest-pressure limit. There is really no reason for it not to be left at its default value of 0. The user may then re-reference to the bottom of the cast, or any intermediate depth, in a separate step.

This is simply copying the behavior of the original Matlab function, which also requires that p_ref be a scalar.

[ledm]

Okay, thanks. I think I may be misunderstanding how this works. Perhaps you can advise me on this issue then. I'm trying to use TEOS-10 to steric calculate sea level rise in CMIP6, based on the work in Paul J Durack et al 2014 Environ. Res. Lett. 9 114017 (https://iopscience.iop.org/article/10.1088/1748-9326/9/11/114017). Although, to be completely honest, I'm struggling with the documentation. While it is not fully described in the paper, they uses EOS-80's function, gpan to calculate dynamic. According to the github page, the TEOS-10 equivalent is geo_strf_dyn_height. However, it's not clear whether they are exactly the same function. Also, the python and R TEOS-10 packages have different descriptions. Are they equivalent? If so, what operator needs to be applied to the TEOS-10 geostrohpic dynicamic height to produce sea surface height?

[ocefpaf]

Closing this one as stale. The defaults are the same in Matlab, Python and R as far and I can tell and a comparison of the EOS-80 gpan and geo_strf_dyn_height is a bit beyond the scope of what the maintainers here can do. @ledm if you figure this out feel free to report back here.