One random thought I wanted to add is this:

I view dask as a pretty advanced python tool. I myself didn't really start to leverage it to its full potential until early last summer I believe. Many students and researchers I work with are still transitioning to xarray and are happy doing things in memory with long runtimes because they view dask as daunting. I think that's totally reasonable -- it was a big jump for me.

So my concern is that when we construct this vectorized grid, if people are working in memory we could end up crashing their server or causing really long runtimes. In some cases maybe looping is better? This goes I guess for the current implementation plus this. But I view @aaronspring's ideas as moving toward a non-looped vectorized system which would create larger objects. I think it'll run beautifully in dask but we can't expect all users to have that.

1. So is there some way to detect their total memory and approximate the size of object we'll generate on the fly?
2. Can we use the looping method when in memory for a certain problem size?

Just some thoughts. We can deal with this in testing.

I see this my way: once people are getting memory issues, they should/need to learn how to use dask to overcome this challenge. if they are way below memory, I dont mind how they do it.

.interp() might be more robust than .sel() if not all times match exactly (offset by 12H) or plain missing dates on obs (although slower). .sel(method='nearest') wouldn't work I think since if model valid time is 2020-04-02, and obs selected the nearest 2020-04-01, it's still a completely different date

So I recommend doing a

try:
    ds.sel()...
except KeyError:
    ds.interp()...
except MemoryError:
    raise MemoryError('Please utilize dask e.g. ds.chunk('auto')')

With interpolation we don’t known which valid time was taken in the end and nan masking stops being effective. So far we required perfect matching of times

in general we should add a multi-dimensional time dimension to our output. quick to calculate with cftime_add_time_from_init_lead this facilitates plotting sooo much:

from climpred.utils import shift_cftime_index,get_lead_cftime_shift_args,shift_cftime_singular
def cftime_add_time_from_init_lead(ds):
    if 'time' not in ds.coords and 'time' not in ds.dims:
        freq=xr.infer_freq(ds.init.to_index())
        time = xr.concat([xr.DataArray(shift_cftime_singular(ds.init,int(lead),freq),dims='init') for lead in ds.lead.values],dim='lead',join='inner',compat='broadcast_equals')    
        time['lead']=ds.lead
        ds.coords['time']=time
    return ds

assert 'time' in cftime_add_time_from_init_lead(fo).coords


hind = cftime_add_time_from_init_lead(hind)

import seaborn as sns
sns.set_palette('viridis',10)
fig,ax=plt.subplots(figsize=(10,6),nrows=2,sharex=False,sharey=True)
for l in mask.lead:
    hind.sel(lead=l)[v].plot(label=str(l.values),x='init',ax=ax[0])
    hind.sel(lead=l)[v].plot(label=str(l.values),x='time',ax=ax[1])
ax[0].set_title('plot by inits')
ax[1].set_title('plot by valid time')
obs_ds_cf_smaller[v].plot(c='k',label='obs',ax=ax[1])
plt.legend()

Oh man, what a great point on plotting.. that was always a pain!

And this will quite easy to add. I will open a PR allowing this soon.

Nevermind; "zero is a zero order spline. It's value at any point is the last raw value seen." not what I thought
https://stackoverflow.com/questions/27698604/what-do-the-different-values-of-the-kind-argument-mean-in-scipy-interpolate-inte

If I am not mistaken, interp(method='zero') is most robust for except clause if you don't want it to be interpolated:

import numpy as np
import pandas as pd
import xarray as xr

cfs_url = 'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/.T/dods'
obs_url = "http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/.monthly/.sst/T/1+index/dods"

obs = xr.open_dataset(obs_url, decode_times=False).isel(T=slice(0, 10))
cfs = xr.open_dataset(cfs_url, decode_times=False)

obs.interp(T=cfs['T'], method='zero')

For missing indices, I think something along the lines of reindex/align; see:
https://xarray.pydata.org/en/stable/indexing.html
but I couldn't get it to work elegantly besides using where, so I'll see if anyone has suggestions:
pydata/xarray#4593

obs = xr.open_dataset(obs_url, decode_times=False).isel(T=slice(0, 10))
cfs = xr.open_dataset(cfs_url, decode_times=False)
tst = obs.interp(T=cfs['T'])
tst.where(tst['T'].isin(obs['T'].values))

I updated the gist and the top post. my new way is faster for non-geospatial but slower for larger than 100000 gridcells.
I think the sel command is tough for big geospatial data. thats really an advantage in the current implementation.

Definitely, we should introduce the new time coord, regardless of what else we do with it.

I'd be interested in your opinions @bradyrx @ahuang11

Lots of code can be used from xarray https://github.com/pydata/xarray/blob/6c32d7c21941461ae9c21b43e6071ee79fb47d68/xarray/core/accessor_dt.py

Could try initialised.coords.to_dataset()