Extracting a single station

[ThomasHaine]

I want to extract and plot a single hydrographic station (i.e., a depth profile at a single latitude and longitude). How do I do this?

E.g., I specify a (lon, lat) pair and plot the model temperature profile at (i) the nearest model grid point or (ii) the specified location using interpolation.

I tried subsample.mooring_array for (i) and subsample.survey_station for (ii) using a single (lon, lat) point. In both cases, it errors, however, saying that more than one point is required.

The docs and tutorial don't show an example of using OceanSpy to plot a single profile, which is a conspicuous omission.

[Mikejmnez]

You can select data locations via the label-based indexes with xarray, including the nearest neighbor. See https://docs.xarray.dev/en/stable/user-guide/indexing.html#indexing

It does not do interpolation. Perhaps oceanspy can provide some help with that, but it might be 2-3 extra lines after selecting the data. Perhaps with should include an example within the documentation on how to do so first, and then asses the whether it is worth it to add it into oceanspy's API

[ThomasHaine]

How would this work? I want to specify the (lon,lat) coordinates, but the fields are distributed over indices like (time, Z, face, Y, X) where X and Y are integer indices. Nearest neighbour model grid point will work fine. Thanks!

[Mikejmnez]

import oceanspy as ospy
import xarray as xr
import numpy as np
from poseidon_viewer import get_shapes
import dask
import xoak

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 12, 8

from dask.distributed import Client
client = Client()

# data from Poseidon viewer 
p1 = [{"type":"Point","coordinates":[-42.04663256856133,52.03406883479701]}]
p2 = [{"type":"Point","coordinates":[-37.736613664969646,52.17340337402146]}]
p3 = [{"type":"Point","coordinates":[-32.2923792604328,49.891501701684035]}]
p4 = [{"type":"Point","coordinates":[-48.625082474043374,48.257101018516096]}]
p5 = [{"type":"Point","coordinates":[-40.00504466686001,47.64941029802745]}]
p6 = [{"type":"Point","coordinates":[-36.148711963646406,46.56872977265587]}]
p7 = [{"type":"Point","coordinates":[-45.222435971207844,44.82612018508269]}]
p8 = [{"type":"Point","coordinates":[-41.592946368183284,41.688286042192715]}]
p9 = [{"type":"Point","coordinates":[-35.92186886345737,42.529707215997405]}]


ps = [p1, p2, p3, p4, p5, p6, p7, p8, p9] # sample as many as you'd like

stations = [ospy.utils.viewer_to_range(ps[i]) for i in range(len(ps))]

# this could be easier if viewer_to_range does it for you....
lons = [stations[i][0][0] for i in range(len(stations))]
lats = [stations[i][1][0] for i in range(len(stations))]

# quick visualization of stations
colors=['r', 'b', 'k', 'orange', 'darkgreen', 'm', 'yellow', 'c', 'indigo']
legends = np.arange(10)
for i in range(len(colors)):
    plt.plot(lons[i], lats[i], 'o', color=colors[i], label='station ' + str(legends[i]))
plt.legend();

# create OceanDataset

od_llc4320 = ospy.open_oceandataset.from_catalog('LLC4320')
od_llc4320._ds = od_llc4320._ds.drop_vars({'k', 'k_u', 'k_p1', 'k_l', 'niter'})
co_list = [var for var in od_llc4320._ds.variables if "time" not in od_llc4320._ds[var].dims]
od_llc4320._ds = od_llc4320._ds.set_coords(co_list)

ds = od_llc4320._ds


varList = ['Temp', 'S', 'U', 'V']
nds = ds.drop_vars([var for var in ds.data_vars if var not in varList])

### here we use xoak. I checked and it doesn't work without it....

# Set lat and lon values are indexers via xoak


if not nds.xoak.index:
    nds.xoak.set_index(["XC", "YC"], index_type = 's2point')


# Create dataset with station data
station = {"XC": ("station", lons), "YC": ("station", lats)}
ds_stns = xr.Dataset(station)


# Select nearest neighbor 

# this takes about 2-3 minutes


stns_dataset = nds.xoak.sel(XC=ds_stns['XC'], YC = ds_stns['YC'])
stns_dataset


# plot vertical profiles...

stns_dataset['Temp'].isel(time=0).squeeze().plot.line(y="Z")
plt.show()

[Mikejmnez]

Note that the example above extracts N different stations from LLC4320, with N (9 in the example above) an arbitrary number.

[malmans2]

It feels like we should support this:

subsample.mooring_array([lon], [lat])

returning a single grid cell with dimension (X, Y, Xp1, Yp1):

              V(x, yp1[1])
U(xp1[0], y)  C(x, y)         U(xp1[1], y)
              V(x, yp1[0])

Then to plot fields other than C, one could use .plot(hue="Xp1"), or select just one face, or taking the average, ...

We already do something much more complicated to get the mooring arrays, getting just one cell should be straightforward (e.g., get the nearest (X,Y) indexes, then infer Xp1 and Yp1 indexes, which is pretty much what is done in the loop to find all moorings).

But I haven't been using OceanSpy in a long time, maybe I'm just overlooking something.

[Mikejmnez]

Thanks @malmans2 ! That is a good point that subsample.mooring_array already does something more complicated and passing single (lat, lon) coordinate should be supported. However oceanspy doesn't support it partly due to some changes I made (in subsample.mooring_array, when computing circle_path trajectory, I require Ymoor to be of size>=2) and some earlier features before my time, for example when computing distance between mooring stations...

I think it is straightforward to allow a single point location to be passed into mooring_array, by adding some extra arguments or conditionals (perhaps something fancier). But in that case we should allow multiple stations to be extracted in parallel and avoid multiple serial evaluations. subsample.mooring_array would connect them via a zig-zag-ing line along a great circle path, not what we want, so we'd have to implement this feature. Not to mention that subsample.mooring_array triggers a cutout either for every single station (if implemented) or the entire subdomain containing the lats/lons passed (this can generate some unnecessary data manipulation when involving LLC-gridded data).

On the other hand, the following code extracts N >=1 arbitrary locations in 4 lines

ds.xoak.set_index(["XC", "YC"], index_type = 's2point')
station = {"XC": ("station", lons), "YC": ("station", lats)}
ds_stns = xr.Dataset(station)
stns_dataset = ds.xoak.sel(XC=ds_stns['XC'], YC = ds_stns['YC'])

the rest of the example code I wrote above is importing, declaring variables and plotting.

I guess my point is that extracting data at a single (or multiple) isolated locations is trivial enough that I don't see much point to implementing it within subsample.mooring_array. Instead, we could just include an example in the documentation or something like it. But that is my opinion and certainly don't discourage any further development on oceanspy.

[ThomasHaine]

How about adding a new function called (something like): subsample.profiles that extracts an ensemble of vertical profiles at the nearest model gridpoints (or interpolates)? To me, that seems most transparent to the user and most aligned with OceanSpy's mission.