something changed in interpolate.cross_section for curvilinear grids

[ahuang11]

Now it's returning all NaNs compared to metpy==0.12.2. Unfortunately I do not have a minimal reproducible example at the moment, but any ideas what might have changed since I don't think the release notes mentioned anything about that being changed? I suspect something to do with metpy.parse_cf() which now is yields a coordinate metpy_crs in my xr.Dataset but still testing.

**curvilinear grid as in 2D lon/lat, 1D x/y

example code

import metpy
from metpy import interpolate
start= [-59.2, 80]
end= [ -42.6, 80]
method = 'nearest'
ds = ds.metpy.parse_cf()
cross = interpolate.cross_section(
    ds, start, end, steps=25, interp_type=method
)

[ahuang11]

I think it has to do with pyproj replacing cartopy.
e4946c6#diff-e0008605917b4b8fb769ea39351ec5161853d1cbe6beb190561b5baa18b555a9

[ahuang11]

0.12.2 yields these x/y

1.0.0 yields these x/y

[jthielen]

Several changes occurred in MetPy 1.0 that could have caused this difference (using PyProj instead of cartopy, and more correctly handling globes/geods). What are the CF parameters for the grid mapping in your file? Based on the values you obtained both before and after, it looks like your points are rather far from your projection origin, which could result in high sensitivities to how the projection is handled. Using the geodesic constructor directly (metpy.interpolate.geodesic) may also help with troubleshooting.

[ahuang11]

The geod differs in versions:
e4946c6#diff-e0008605917b4b8fb769ea39351ec5161853d1cbe6beb190561b5baa18b555a9R101-R106
print(g, p)
print(crs.proj4_init)

1.0.0

Geod('+a=6378137 +f=0.0033528106647475126')
proj=lcc lat_0=0 lon_0=70.127 lat_1=-30 lat_2=-60 x_0=0 y_0=0 datum=WGS84 units=m no_defs ellps=WGS84 towgs84=0,0,0

0.12.2

Geod(ellps='sphere')
+ellps=sphere +proj=lcc +lon_0=70.127 +lat_0=39.0 +x_0=0.0 +y_0=0.0 +lat_1=-30 +lat_2=-60 +no_defs

[ahuang11]

I don't think the Geod is initialized properly. It's missing the proper standard parallels and central lons.

grid_mapping_name :
lambert_conformal_conic
standard_parallel :
[-30, -60]
central_longitude :
70.1273
longitude_of_central_meridian :
70.127
long_name :
coordinate reference system

[jthielen]

Hmm...so in the 1.0 verison, what is pyproj.Proj(ds.metpy.pyproj_crs)?

[ahuang11]

Description: PROJ-based coordinate operation Area of Use: - undefined -

[ahuang11]

Maybe my attrs are defined incorrectly compared to the example NARR:

grid_mapping_name :
lambert_conformal_conic
standard_parallel :
50.0
longitude_of_central_meridian :
-107.0
latitude_of_projection_origin :
50.0

Mine:

grid_mapping_name :
lambert_conformal_conic
standard_parallel :
[-30, -60]
central_longitude :
70.1273
longitude_of_central_meridian :
70.127
long_name :
coordinate reference system

Nevermind, it looks like it's properly named...

[jthielen]

It does look properly named...but now that you point it out, those values with the example NARR look suspiciously like the ones in #1646 (comment). Would you be able to verify that your script is not mistakenly associating the wrong dataset or wrong dataset's CRS?

[ahuang11]

You're right, my bad; I am jumping around in my notebook; I updated my comment and lat_0 seems to be off.

Geod('+a=6378137 +f=0.0033528106647475126')
proj=lcc lat_0=0 lon_0=70.127 lat_1=-30 lat_2=-60 x_0=0 y_0=0 datum=WGS84 units=m no_defs ellps=WGS84 towgs84=0,0,0

[ahuang11]

Okay I manually added
ds['crs'].attrs['latitude_of_projection_origin'] = 39.0
But still getting NaNs.

Geod('+a=6378137 +f=0.0033528106647475126')
proj=lcc lat_0=39 lon_0=70.127 lat_1=-30 lat_2=-60 x_0=0 y_0=0 datum=WGS84 units=m no_defs ellps=WGS84 towgs84=0,0,0

Is there a difference betwteen
ellps=WGS84
and
ellps=sphere

The only differences I see now are:

Geod(ellps='sphere')
ellps='sphere'

Geod('+a=6378137 +f=0.0033528106647475126')
datum=WGS84 towgs84=0,0,0 ellps=WGS84

[jthielen]

Yes, there is indeed a difference between the two, as WGS84 non-spherical (it has a flattening, meaning major and minor axes are different). How different are your x and y coordinates now (or the geodesic path constructed by metpy.interpolate.geodesic) after the addition of the origin latitude argument?

[ahuang11]

I added a few more attrs, but still off

ds['crs'].attrs['latitude_of_projection_origin'] = 39.0
ds['crs'].attrs['earth_shape'] = 'spherical'
ds['crs'].attrs['earth_radius'] = 6367470.21484375

1.0.0

Geod('+a=6367470.21484375 +f=0')
proj=lcc lat_0=39 lon_0=70.127 lat_1=-30 lat_2=-60 x_0=0 y_0=0 R=6367470.21484375 units=m no_defs

0.12.2

[ahuang11]

Is there a way to force
ellps='sphere'

[jthielen]

I think the default radius for ellps='sphere' is 6370997, so would you be able to try that instead to replicate it? I'm not sure there is a way for force ellps='sphere' through pyproj.CRS.from_cf (which is all MetPy's using under-the-hood)

[ahuang11]

Oh it matches now.

Ah off bounds. I was testing the wrong latitudes/longs

[ahuang11]

Awesome thanks for helping me debug this! Apparently it's just my attributes too dependent on cartopy's defaults!

[ahuang11]

Hmm, if I set latitude_of_projection_origin, I get

da['metpy_crs'].item().to_cartopy()

    return ccrs.Mercator(globe=globe, **kwargs)
TypeError: __init__() got an unexpected keyword argument 'central_latitude'

[ahuang11]

Eh maybe I just have to pop the added kwarg if == mercator