Reference pressure array in geostrophy.geo_strf_dyn_height
Closed this issue · 7 comments
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.
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
"""
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.
prefv = p_ref[ind[1]]I think you are making an assumption here about the dimensions of the input and the value of 'axis'.
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.
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.
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?