Can not execute cross_section - 2D LATS and LONS are in the separate file

[VanuatuN]

Dear MetPy team,

I've spent several hours today trying to plot cross section from NEMO ocean model output, but there are still issues.
I would be very grateful if you can advice me how to proceed with it.

My 2D lon and lats are stored in a separate file, like this:

My 3D data is in another file, like this:

The second file also contain nav_lat and nav_lon, but there is an issue with missing values in coors, so I need to use lats
and lons from the separate file.

I've tried to execute https://unidata.github.io/MetPy/latest/examples/cross_section.html#sphx-glr-examples-cross-section-py

example, but obviously got error message.

The thing is that I only have 2D lats and lons and no projection naming is available netcdf desscription.
This is how it looks:

Thank you very much for any advise!!!

Warm regards,
Natalia

[jthielen]

Given the current way that xarray does coordinate indexing, there are two options to explore for taking a cross section when you just have 2D longitude and latitude and not dimension coordinates from your source data:

1) Add in the dimension coordinates

By using .metpy.assign_y_x(), this is what it looks like you were trying to do, so you were on the right track! However, the dataset at that point doesn't have the metpy_crs coordinate (as noted in the error message), so there isn't enough for that method to do its job. Given that it doesn't sound like the grid projection metadata is present in these files (so, can't use parse_cf()), you would need to obtain this from elsewhere (and so use assign_crs()). From the plots, these look Mercator-like, so that's what I'll demo below, but you should update them to the proper values (following CF).

# Assign the CRS info
data = data.metpy.assign_crs({
    'grid_mapping_name': 'mercator',
    'longitude_of_projection_origin': 0,
    'standard_parallel': 0
}

# Include dimension coordinates
data = data.metpy.assign_y_x()

# Proceed with cross section
start = (30, -77)
end = (45, -64)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))

2) Use xoak to directly, but slowly, index with longitude and latitude

If you have no ability to get the true CRS information and are therefore stuck with just 2D longitude and latitude, then your next best option is to use xoak, which supports direct indexing of multi-dimensional longitude and latitude. This uses a search tree, so it's naturally slower than 1D indexing with dimension coordinates, but it's the best alternative when you don't have suitable 1D indexes (and much, much better than Python-level looping to search through the coordinates). An example can be found here: https://pop-tools.readthedocs.io/en/latest/examples/xoak-example.html. Otherwise, to replicate what MetPy's cross section does (albeit with nearest-neighbor interpolation, which is the only available option, rather than MetPy's default linear interpolation), it could look something like

from metpy.interpolate import geodesic
import pyproj
import xoak

# Set up search index
data.xoak.set_index(['nav_lat', 'nav_lon'], 'scipy_kdtree')

# Set up path of the cross section
wgs84 = pyproj.CRS(4326)
start = (30, -77)
end = (45, -64)
path = geodesic(wgs84, start, end, steps=100)

# Take the cross section
cross = data.xoak.sel(
    nav_lat=xr.DataArray(cross[:, 1], dims='index', attrs=data.nav_lat.attrs),
    nav_lon=xr.DataArray(cross[:, 0], dims='index', attrs=data.nav_lon.attrs)
)

(note: the default value of steps at 100 will likely cause artifacts due to the nearest neighbor interpolation. You'll likely need to adjust this to a better value given your path length and data resolution.)

[VanuatuN]

Dear jthielen,

Thank you very much for the answer!

I was trying to did exactly what you are suggested as the first option and got the following:

Would be very extremely grateful for any suggestions....

Thank you and warm regards,
Natalia

[VanuatuN]

Ah, sorry, I haven't updated characteristics of the projection..
Now need to sort out how to find these
'longitude_of_projection_origin': 0,
'standard_parallel': 0

Will try to investigate into this..

[VanuatuN]

Guys, I did it, so happy! Here is the code that worked in my case, many thanks for help!!!

from metpy.interpolate import geodesic
import pyproj
import xoak

### Set up search index
data.xoak.set_index(['nav_lat', 'nav_lon'], 'scipy_kdtree')

### Set up path of the cross section
wgs84 = pyproj.CRS(4326)
start = (30, -77)
end = (45, -64)
path = geodesic(wgs84, start, end, steps=100)
path
   
### Take the cross section
cros = data.xoak.sel(
                   nav_lat=xr.DataArray(path[:, 1], dims='index', attrs=coors.nav_lat.attrs),
                   nav_lon=xr.DataArray(path[:, 0], dims='index', attrs=coors.nav_lon.attrs)
                    ) 
cros