def cars_profile(filename, doy, latitude, longitude, depth):
"""
For now only the nearest value
For now only for one position, not an array of positions
longitude 0-360
"""
assert np.size(doy) == 1
assert np.size(latitude) == 1
assert np.size(longitude) == 1
#assert np.size(depth) == 1
assert (longitude >= 0) & (longitude <= 360)
assert depth >= 0
nc = netCDF4.Dataset(filename)
t = 2 * np.pi * doy/366
# Improve this. Optimize to get only necessary Z
Z = slice(0, nc['depth'].size)
I = np.absolute(nc['lat'][:] - latitude).argmin()
J = np.absolute(nc['lon'][:] - longitude).argmin()
# Not efficient, but it works
assert (nc['depth'][:64] == nc['depth_ann'][:]).all()
assert (nc['depth'][:55] == nc['depth_semiann'][:]).all()
value = nc['mean'][:, I, J]
value[:64] += nc['an_cos'][Z, I, J] * np.cos(t) + \
nc['an_sin'][:, I, J] * np.sin(t)
value[:55] += nc['sa_cos'][Z, I, J] * np.cos(2*t) + \
nc['sa_sin'][:, I, J] * np.sin(2*t)
value = value
output = {'depth': np.asanyarray(depth)}
from scipy.interpolate import griddata
output['value'] = griddata(nc['depth'][Z], value[Z], depth)
for v in ['std_dev']:
output[v] = griddata(nc['depth'][Z], nc[v][Z, I, J], depth)
return output
评论列表
文章目录