Taking a cross section of data (2d lat and lon)

[pecos27]

I am using this model data output for salinity and temperature, and it has 2d lons and lats, and I am hoping to take a cross section of the data. Basically to get an output array of the field of interest, "ds_along_section", which has the "nearest neighbors" along the cross section lined up along a new dimension. The new dimension "ds_along_section" later represents distance which can be calculated using the Haversine formula. The data array has a shape of time, depth, y, and x but has coordinates lon and lat which are both 2-dimensional.

For 1D lon and latitude, xarray's vectorized selections works wonders, but with 2D lat and lon arrays there becomes a limitation.
Does anyone have any recommendations on selecting a cross section over 2D lat and lon?
I tried the package MetPy which performs a cross sectional analysis, but the data does not appear to fit MetPy's data structure parameters.

Unfortunately, this is not so reproducible, as the model data is needed. As an additional attempt to take a cross section, I took a cross section using the package region mask:

ds = xr.open_dataset(inpath+'infile')

small_area = regionmask.Regions([
        np.array([[765., 1110], [766, 1111], [826, 1021],[825, 1020]]) 
    ],
    names=["area", ])

cross_sec = salinity_1993_2019.where(~small_area.mask(salinity_1993_2019).astype(bool),drop=True)

cross_sec.shape

(12, 46, 90, 61) #time, depth, y, and x (y being the indices from 0, 1,... end, likewise with x)

cross_sec[0,0].plot(); #[0,0] being first time index and then the surface layer 

Is it possible to drop all nans except along the "line" of data? If that makes sense?
I tried dropna for both xarray data arrays and pandas data frames, but then the entire axis was dropped (not just the nan values).
This is just an attempt at taking a diagonal slice across the temp/salinity... would this be a viable way to take a cross section?
Or is there a more efficient/better way to take a cross section over 2D lat and lon arrays?
Thank you! Hopefully the question is understandable (any tips or ideas for taking a cross section would be greatly appreciated). Apologies for the lack of reproducibility.

[jthielen]

As you noted, having 1D dimension coordinates enables xarray's vectorized indexing functionality (which is all MetPy is doing under the hood). And so, one solution to taking a cross section when you have 2D lon and lat arrays (and probably the most efficient) is to re-derive the "missing" x and y 1D dimension coordinates given your data's projection (if sticking with MetPy, assign_y_x should do the trick if you've added the projection information properly), and then use that built-in vectorized indexing.

Alternatively, you can index the data with the 2D lon and lat arrays directly using the extension package xoak, as in this example. While this is likely to be less efficient (and may require using Dask if your dataset is too large), it may be the only option if you don't have sufficient information on your data projection or you have a complex grid that isn't handled by pyproj or cartopy.