cross_section() does not match requested start and end

[sgdecker]

Here is my code:

import pickle
import numpy as np
import metpy
from metpy.interpolate import cross_section
import xarray as xr
import cartopy
import cartopy.crs as ccrs

with open('ds2.pickle', 'rb') as f:
    ds2 = pickle.load(f)

start = (40.0, -105.0)
end = (45.0, -68.0)

cross = cross_section(ds2, start, end)
cross.set_coords(('lat', 'lon'), True)

print(cross.lat[0].values, cross.lon[0].values)
print(cross.lat[-1].values, cross.lon[-1].values)

# Extra printouts to verify projection information agreed
#print(ds2.crs.to_dict()['data'].to_dict())
#print(ds2.LambertConformal)

print(metpy.__version__)
print(xr.__version__)
print(cartopy.__version__)
 

The output:

xsecbug.py:16: FutureWarning: The inplace argument has been deprecated and will be removed in a future version of xarray.
  cross.set_coords(('lat', 'lon'), True)
42.14027577619222 -91.72784705269194
43.045485539202936 -54.08897522118793
0.10.2
0.12.2
0.17.0

• Expected output
  I would hope the start and end points in the cross section would be much closer to the (40, -105) and (45,-68) requested.

The pickle file can be downloaded from https://drive.google.com/file/d/1B-e1FIZMh6T6Px40FCtcfiGzKVxy8Zof/view?usp=sharing

[jthielen]

Would you be able to determine how the original lat/lon values in dataset were computed (and on what datum)?

As far as I can tell looking at the dataset, they don't seem to match the x/y dimension coordinates and crs given (which are what MetPy uses for finding the cross section) under cartopy transforms. For example, it looks like the longitude values given are off by around 10-17 degrees from what would be expected from the x/y coordinates and crs:

crs = cross['geopotential_height'].metpy.cartopy_crs
globe = cross['geopotential_height'].metpy.cartopy_globe

x, y = xr.broadcast(ds2.x, ds2.y)
lonlats = ccrs.Geodetic(globe).transform_points(crs, x.values, y.values)
lons = lonlats[..., 0].transpose()
lats = lonlats[..., 1].transpose()

(ds2.lon - lons).plot()

As a sanity check, the computed lat/lons from the cross section match the desired start and end points:

for i in (0, -1):
    lon, lat = ccrs.Geodetic(globe).transform_point(cross.x[i], cross.y[i], crs)
    print(lat, lon)

40.00000000000001 -105.00000000000001
45.0 -68.0

[sgdecker]

OK, I figured something wonky was going on in the coordinate calculations. This is WRF output that has been post-processed. What's weird is a different WRF output file (albeit for a different domain), when processed in the same way, generates the proper cross section. I'll need to spend some time reducing the post-processing piece to its essence (right now the "bad" WRF output is a 44 GB netCDF file), and then I'll post it here.

In the meantime, I am wondering if it might be a good idea for cross_section to perform some of these sanity checks and error out when confronted with a DataArray/Dataset like this.

[sgdecker]

Here's an updated testcase:

import numpy as np
import netCDF4 as nc
import metpy
import metpy.calc as mpcalc
from metpy.units import units
from metpy.interpolate import log_interpolate_1d, cross_section
import xarray as xr
import pyproj
import cartopy

def av_ar(a, n=1, axis=-1):
	a = np.asanyarray(a)
	nd = a.ndim
	axis = np.core.multiarray.normalize_axis_index(axis, nd)
	slice1 = [slice(None)] * nd
	slice2 = [slice(None)] * nd
	slice1[axis] = slice(1, None)
	slice2[axis] = slice(None, -1)
	slice1 = tuple(slice1)
	slice2 = tuple(slice2)
	op = not_equal if a.dtype == np.bool_ else np.add
	for _ in range(n):
		a = op(a[slice1], a[slice2]) / 2.
	return a

def read_var(NCfile, varname, hr):
	fid = nc.Dataset(NCfile, 'r')
	var_out = fid.variables[varname][:]
	var_time = var_out[hr]
	fid.close()
	return var_time

def read_map_parm(NCfile):
	fid = nc.Dataset(NCfile, 'r')
	dx_val, dy_val = fid.DX, fid.DY
	nx_val, ny_val = fid.dimensions['west_east'].size, fid.dimensions['south_north'].size
	y1, y2, y0, ym, x0 = fid.TRUELAT1, fid.TRUELAT2, fid.MOAD_CEN_LAT, fid.CEN_LAT, fid.STAND_LON
	fid.close()
	return y1, y2, y0, ym, x0, dy_val, dx_val, ny_val, nx_val

def calc_pres(fid_name, hr):
	p0 = read_var(fid_name, 'PB', hr) * units.Pa
	p1 = read_var(fid_name, 'P', hr) * units.Pa
	pres_pa = p0 + p1
	return pres_pa.to('hPa')

def calc_hght(fid_name, hr):
	gvty0 = read_var(fid_name, 'PHB', hr)
	gvty1 = read_var(fid_name, 'PH', hr)
	geopot0 = gvty0 + gvty1
	geopot1 = av_ar(geopot0, axis=0) * units.m**2 / units.s**2
	h = mpcalc.geopotential_to_height(geopot1)
	return h

def create_deltas(fid_name, dx_in):
        # assume same grid spacing for both x and y directions
        # make sure array size is 1 less than data in dir of interest
	map_x = read_var(fid_name, 'MAPFAC_U', 0)[:,1:-1]
	map_y = read_var(fid_name, 'MAPFAC_V', 0)[1:-1,:]
	dx_out = np.empty_like(map_x)
	dy_out = np.empty_like(map_y)
	dx_out = dx_in / map_x
	dy_out = dx_in / map_y
	return dy_out, dx_out # following same order convention as wrf output

def sigma2iso(p, h, th, dx, dy, l):
	plevs = np.arange(1025., 40., -25.) * units.hPa
	h_int, thta_int = log_interpolate_1d(plevs, p, h, th, axis=0)
	return plevs, h_int, thta_int

f_input = 'wrfout1.nc'
h_ind = 1

yp1, yp2, yp0, yp0_m, xp0, dy, dx, ny, nx = read_map_parm(f_input)
lats = read_var(f_input, 'XLAT', 0)
lons = read_var(f_input, 'XLONG', 0)
thta = (read_var(f_input, 'T', h_ind) + 300.0) * units.K
pres = calc_pres(f_input, h_ind)
hght = calc_hght(f_input, h_ind)
dy_var, dx_var = create_deltas(f_input, dy*units.m)
plevs, h_iso, th_iso = sigma2iso(pres, hght, thta, dx_var, dy_var, lats)

# Define the WRF projection
wrf_proj = pyproj.Proj(proj='lcc',
                       lat_1=yp1, lat_2=yp2,
                       lat_0=yp0, lon_0=xp0,
                       a=6370000, b=6370000)

# Easting and Northing of the domain center point
wgs_proj = pyproj.Proj(proj='latlong', datum='WGS84')
e, n = pyproj.transform(wgs_proj, wrf_proj, xp0, yp0_m)

# Lower left corner of the domain
y0 = -(ny-1) / 2. * dy + n
x0 = -(nx-1) / 2. * dx + e

# Get grid values
x, y = np.arange(nx) * dx + x0, np.arange(ny) * dy + y0

# create the xarray Dataset
ds2 = xr.Dataset({'potential_temperature': (['vertical', 'y', 'x'], th_iso,
                                            {'units': str(th_iso.units)}),
				'geopotential_height': (['vertical', 'y', 'x'], h_iso,
                                                        {'units': str(h_iso.units)}),
				'lon': (['y', 'x'], lons, {'units': 'degrees_east'}),
				'lat': (['y', 'x'], lats, {'units': 'degrees_north'})},
				coords={'vertical': (['vertical'], plevs,
                                                     {'units': str(plevs.units)}),
				'y': (['y'], y, {'units': 'meter'}),
				'x': (['x'], x, {'units': 'meter'})})

# Define the grid mapping
ds2['LambertConformal'] = xr.DataArray(np.array(0), attrs={
	'grid_mapping_name': 'lambert_conformal_conic',
	'earth_radius': 6370000,
	'standard_parallel': (yp1, yp2),
	'longitude_of_central_meridian': xp0,
	'latitude_of_projection_origin': yp0})

for var in ds2.data_vars:
	ds2[var].attrs['grid_mapping'] = 'LambertConformal'
	if 'projection' in ds2[var].attrs:
		del ds2[var].attrs['projection']
	if 'coordinates' in ds2[var].attrs:
		del ds2[var].attrs['coordinates']

ds2 = ds2.metpy.parse_cf().squeeze()

start = (40.0, -105.0)
end = (45.0, -82.0)

cross = cross_section(ds2, start, end)
cross.set_coords(('lat', 'lon'), True)

print(cross.lat[0].values, cross.lon[0].values)
print(cross.lat[-1].values, cross.lon[-1].values)

print(metpy.__version__)
print(xr.__version__)
print(cartopy.__version__)

wrfout1.nc fails, but change line 71 to refer to wrfout2.nc and it works as expected. The two netCDF files are at
https://drive.google.com/file/d/1VbpkIKsEX2tHHi07KCeDuJ1iFq3G6xSJ/view?usp=sharing
and
https://drive.google.com/file/d/1Vd5HLck8R3Md11gVSiOBsLsHflE2-4Cs/view?usp=sharing
respectively.

[jthielen]

I definitely like the idea of being able to sanity check the coordinates of a DataArray/Dataset, but I hesitate in making it tied to cross_section specifically, since any coordinate errors will not just affect cross_section, but also the kinematics functions in the near future (ref #893). Unless there are suggestions to the contrary, I think I'll just try including something like this as I work on dataset helpers.

This would also require identification of both lat/lon auxiliary coordinates and x/y dimension coordinates at the same time, which the current implementation of coordinate identification cannot do. But, I think it would be a worthwhile change (see #1090).

[sgdecker]

Aha, I see the WRF output that works has CEN_LON and STAND_LON attributes that are the same, but for the wonky output, CEN_LON is -72.26288, but STAND_LON is -85.5. This is suspiciously within the range of 10-17 degrees found by @jthielen above. I will cross my fingers that x0 in my posted code should be CEN_LON, not STAND_LON, but I suppose it is possible I need CEN_LON for some parts of the calculation, and STAND_LON for others. I'll find out soon enough!

[sgdecker]

It's the latter case (and I should have written xp0 above). Looks like I need both STAND_LON and CEN_LON, and I'll need to figure out where to use each one.

[sgdecker]

OK, I think I have it figured out. The relevant sections of the code should be modified as follows, where xs is STAND_LON, and xp0 is CEN_LON:

# The pyproj object should use the STAND_LON value:
wrf_proj = pyproj.Proj(proj='lcc',
                       lat_1=yp1, lat_2=yp2,
                       lat_0=yp0, lon_0=xs,
                       a=6370000, b=6370000)

To compute easting and northing, use CEN_LON:

e, n = pyproj.transform(wgs_proj, wrf_proj, xp0, yp0_m)

The xarray metadata should use STAND_LON (consistent with pyproj):

ds2['LambertConformal'] = xr.DataArray(np.array(0), attrs={
	'grid_mapping_name': 'lambert_conformal_conic',
	'earth_radius': 6370000,
	'standard_parallel': (yp1, yp2),
	'longitude_of_central_meridian': xs,
	'latitude_of_projection_origin': yp0})

[dopplershift]

Looks like things have been figured out, but we'll leave this open as a reminder for more sanity checking in the cross section code.