Skip to content

cross_section() does not match requested start and end #1089

New issue
New issue
Open
Open
cross_section() does not match requested start and end#1089
Labels
Area: Cross-sectionsPertains to making cross-sections through dataType: EnhancementEnhancement to existing functionality

Description

@sgdecker
sgdecker
opened on Jul 10, 2019

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

Activity

  1. jthielen commented on Jul 10, 2019

    @jthielen
    jthielen
    on Jul 10, 2019
    Collaborator

    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

  2. sgdecker commented on Jul 10, 2019

    @sgdecker
    sgdecker
    on Jul 10, 2019
    ContributorAuthor

    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.

  3. sgdecker commented on Jul 10, 2019

    @sgdecker
    sgdecker
    on Jul 10, 2019
    ContributorAuthor

    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.

  4. jthielen commented on Jul 10, 2019

    @jthielen
    jthielen
    on Jul 10, 2019
    Collaborator

    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).

  5. sgdecker commented on Jul 11, 2019

    @sgdecker
    sgdecker
    on Jul 11, 2019
    ContributorAuthor

    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!

  6. sgdecker commented on Jul 11, 2019

    @sgdecker
    sgdecker
    on Jul 11, 2019
    ContributorAuthor

    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.

  7. sgdecker commented on Jul 11, 2019

    @sgdecker
    sgdecker
    on Jul 11, 2019
    ContributorAuthor

    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})

  8. dopplershift commented on Jul 26, 2019

    @dopplershift
    dopplershift
    on Jul 26, 2019
    Member

    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.

  9. dopplershift
    added
    Area: Cross-sectionsPertains to making cross-sections through data
    Type: EnhancementEnhancement to existing functionality
    on Jul 26, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    Area: Cross-sectionsPertains to making cross-sections through dataType: EnhancementEnhancement to existing functionality

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions