hainegroup / oceanspy

A Python package to facilitate ocean model data analysis and visualization.
https://oceanspy.readthedocs.io
MIT License
96 stars 32 forks source link

Order of stations error in `subsample_stations` #375

Closed ThomasHaine closed 1 year ago

ThomasHaine commented 1 year ago

@Mikejmnez the subsample_stations method returns stations ordered backwards. Try:

import oceanspy as ospy
from poseidon_viewer import get_shapes
import numpy as np
import xoak
import matplotlib.pyplot as plt

od = ospy.open_oceandataset.from_catalog('ECCO')

# From Poseidon-viewer:
ps = [{"type":"Point","coordinates":[-56.761672775466664,60.13158802096402]},{"type":"Point","coordinates":[-39.43864698533882,23.746791019792084]}]

# Create dataset with station data for the time of interest
lons, lats = ospy.utils.viewer_to_range(ps)
stations_ds = od.subsample.stations(Ycoords=lats, Xcoords=lons)
stations_ds._ds = stations_ds._ds.sel(time = '1992-03-16T12:00:00.000000000')
Nstations = len(lons)

# Plot station locations...
cut_od = od.subsample.cutout(XRange = [min(lons)-10,max(lons)+10],YRange = [min(lats)-5,max(lats)+5])
fig = plt.figure(figsize=(12, 5))
ax = cut_od.plot.horizontal_section(varName="Depth", cmap='Greys_r')
colors=['r', 'b', 'k', 'orange', 'darkgreen', 'm', 'yellow', 'c', 'indigo', 'magenta', 'green']
legestations_ds = np.arange(Nstations)
for i in range(Nstations):
    plt.plot(lons[i], lats[i], 'o', color=colors[i], label='station ' + str(legestations_ds[i]))
plt.legend();
plt.show()

# plot vertical profiles...
for i in range(Nstations):

    plt.subplot(1,2,1)
    plt.plot(stations_ds["THETA"].isel(station = i),stations_ds["Z"], color=colors[i], label='station ' + str(legestations_ds[i]))

    plt.subplot(1,2,2)
    plt.plot(stations_ds["SALT"].isel(station = i),stations_ds["Z"], color=colors[i], label='station ' + str(legestations_ds[i]))

plt.legend();
plt.show()

One station (0, red) is in the Labrador Sea, the other in the subtropical Atlantic (1, blue). But the vertical profiles have them backwards (see screenshot and look at temperature profiles on the left).

image
Mikejmnez commented 1 year ago

Important info about this sampleMethod:

The code you are referring in this issue is only accessible through a branch on my local repo (which you forked and pip installed), meaning this sampleMethod is not yet accessible via hainegroup/oceanspy. This particular branch and therefore sampleMethod is under active development, and broadly remains untested... I am getting close to finishing but there are details that remain unfixed (like accurate representation of velocity/corner points in the arctic)

With that said I tried to reproduce the behavior you describe with my up-to-date fork-branch and couldn't...

%%time

# From Poseidon-viewer:
ps = [{"type":"Point","coordinates":[-56.761672775466664, 60.13158802096402]},{"type":"Point","coordinates":[-39.43864698533882,23.746791019792084]}]

# Create dataset with station data for the time of interest
lons, lats = ospy.utils.viewer_to_range(ps)
Nstations = len(lons)

args = {
    "Xcoords": lons, 
    "Ycoords": lats,
    "tcoords":'1992-03-16T12:00:00.000000000'
}
od_stations = od.subsample.stations(**args)
cut_od = od.subsample.cutout(XRange = [min(lons)-10,max(lons)+10],YRange = [min(lats)-5,max(lats)+5])
fig = plt.figure(figsize=(12, 5))
ax = cut_od.plot.horizontal_section(varName="Depth", cmap='Greys_r')
colors=['r', 'b', 'k', 'orange', 'darkgreen', 'm', 'yellow', 'c', 'indigo', 'magenta', 'green']
legestations_ds = np.arange(Nstations)
for i in range(Nstations):
    plt.plot(lons[i], lats[i], 'o', color=colors[i], label='station ' + str(legestations_ds[i]))
plt.legend();
plt.show()
# plot vertical profiles...
for i in range(Nstations):
    plt.subplot(1,2,1)
    plt.plot(od_stations._ds["THETA"].isel(station = i).squeeze(), od_stations._ds["Z"], color=colors[i], label='station ' + str(legestations_ds[i]))

    plt.subplot(1,2,2)
    plt.plot(od_stations._ds["SALT"].isel(station = i).squeeze(), od_stations._ds["Z"], color=colors[i], label='station ' + str(legestations_ds[i]))
plt.legend();
plt.show()

stations

ThomasHaine commented 1 year ago

Thanks @Mikejmnez. You're right, I should have raised this on your fork. Updating to your latest version fixes the issue!