`survey_stations` station locations padded with zeros for `ECCO`

[ThomasHaine]

I'm having an issue with the ECCO dataset, possibly related to the landmask. If I make a cutout then use survey_stations with LLC4320 or EGshelfIIseas2km_ASR_full and then inspect the XC, YC coordinates, all makes sense. With ECCO, I get trailing zeros on XC, YC, however. See the SciServer MWE below. What's going on?

print(ospy.__version__)      # Version information
from dask.distributed import Client
client = Client()

OGCM_dataset = 'ECCO'
# OGCM_dataset = 'LLC4320'
# OGCM_dataset = 'EGshelfIIseas2km_ASR_full'
od = ospy.open_oceandataset.from_catalog(OGCM_dataset)

if OGCM_dataset == 'LLC4320':
    od._ds = od._ds.drop_vars({'k', 'k_u', 'k_p1', 'k_l'})
    co_list = [var for var in od._ds.variables if "time" not in od._ds[var].dims]
    od._ds = od._ds.set_coords(co_list)
elif OGCM_dataset == 'ECCO':
    od._ds = od._ds.drop_vars('time_bnds')
    ds = od.dataset
    ds['Temp'] = ds['THETA']
    od = ospy.OceanDataset(ds)

varList = ['Temp']        # Specify variables of interest here.

this_time = '2012-04-25T00'
lats = [   60.0699,
   60.0477,
   60.0307,
   60.0061,
   59.9880,
   59.9579,
   59.9265,
   59.9041,
   59.8610,
   59.8158,
   59.5808,
   59.2147,
   59.0988,
   59.0205,
   58.9555,
   58.8853]
lons = [  -42.8257,
  -42.5995,
  -42.4285,
  -42.2053,
  -42.0259,
  -41.7395,
  -41.4377,
  -41.1401,
  -40.6870,
  -40.2786,
  -37.7983,
  -35.1258,
  -33.6822,
  -32.7675,
  -31.9590,
  -31.3030]

cutout_kwargs = {
    'varList': varList,
    'timeRange': this_time,
    'XRange': lons, 
    'YRange': lats, 
    'add_Hbdr': True,
}

# Cutout region of interest:
cut_od = od.subsample.cutout(**cutout_kwargs)

# Extract survey stations
od_surv = cut_od.subsample.survey_stations(
    Xsurv=lons,
    Ysurv=lats,
    ZRange=[0,-4000],  # controls deepest depth on plots
    delta=10.0  # survey station spacing in km
)
od_surv.dataset['YC'].values

[Mikejmnez]

I will take a look at this during the week.

[Mikejmnez]

I finally had some time to figure out a way to get around this. The issue is likely due to the coarse resolution of ECCO and the interpolation near the edge of the dataset after performing a cutout.

workaround this.

In ECCO and similar LLC datasets, there is no edge/boundary, and so this shouldn't be an issue. However, the reason it's happening is because you are doing a double cutout to get the survey station. By default, when performing a survey_station (or mooring_array) subsampling there is an internal call to cutout. During this internal call, XRange and YRange are byX_surv and Y_surv, along with a narrow horizontal boundary Hbdr increase.

Because of default behavior, there is no need to first manually perform acutout and then survey_station (or mooring_array if that were the case), since oceanspy is already doing that for you. When you do a cutout first, you reduce the domain to the size of the survey track, and when you then compute the survey_station, the coarse interpolation then interacts with the edge of the dataset (I still need to spend some time to better understand this), which results in nans or zeros at the right edge of the track.

workaround: Survey_station from original OceanDataset.

Keeping everything in your example the same except that when computing the survey you do so from the original od.

# Extract survey stations
od_surv = od.subsample.survey_stations(
    Xsurv=lons,
    Ysurv=lats,
    varList=varList,
    ZRange=[0, -4000],  # controls deepest depth on plots
    delta=20  # survey station spacing in km
)

Then you plot

ax1 = cut_od.plot.horizontal_section(varName="YC", cmap='Greys_r')
Xm, Ym = od_surv.dataset['XC'].values, od_surv.dataset['YC'].values
ax1.plot(Xm, Ym, 'r', lw=3, marker='d', markersize=5);

you get

[ThomasHaine]

Great! All works as expected now.