nilyin
/
wrf-python
forked from 3rdparty/wrf-python
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Fixes vertcross and interpline when using lat/lon coordinates with a projection object. 

The fix requires the user to supply the lower left lat/lon corner points for the domain. Added unit tests. Updated docs. Fixes #44.
lon0
Bill Ladwig 7 years ago
parent
f02db010f8
commit
b8bc1fbb38
5 changed files with 96 additions and 20 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 12
      doc/source/basic_usage.rst
  2. 30
      src/wrf/interp.py
  3. 24
      src/wrf/interputils.py
  4. 18
      src/wrf/metadecorators.py
  5. 32
      test/utests.py

12
doc/source/basic_usage.rst
Unescape View File

@ -736,8 +736,9 @@ Example Using Lat/Lon Coordinates
start_point = CoordPair(lat=start_lat, lon=start_lon) start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon) end_point = CoordPair(lat=end_lat, lon=end_lon)
# When using lat/lon coordinates, you must supply a netcdf file object, or a # When using lat/lon coordinates, you must supply a WRF netcdf file object,
# projection object. # or a projection object with the lower left latitude and lower left
# longitude points.
p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True) p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)
print(p_vert) print(p_vert)
@ -983,8 +984,11 @@ Example Using Lat/Lon Coordinates
start_point = CoordPair(lat=start_lat, lon=start_lon) start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon) end_point = CoordPair(lat=end_lat, lon=end_lon)
# Calculate the vertical cross section. By setting latlon to True, this # Calculate the interpolated line. To use latitude and longitude points,
# also calculates the latitude and longitude coordinates along the line # you must supply a WRF NetCDF file object, or a projection object along
# with the lower left latitude and lower left longitude points.
# Also, by setting latlon to True, this routine will
# also calculate the latitude and longitude coordinates along the line
# and adds them to the metadata to help with plotting labels. # and adds them to the metadata to help with plotting labels.
t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True) t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)

30
src/wrf/interp.py
Unescape View File

@ -91,7 +91,8 @@ def interplevel(field3d, vert, desiredlev, missing=default_fill(np.float64),
@set_interp_metadata("cross") @set_interp_metadata("cross")
def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
wrfin=None, timeidx=0, stagger=None, projection=None, wrfin=None, timeidx=0, stagger=None, projection=None,
ll_lat=None, ll_lon=None,
pivot_point=None, angle=None, pivot_point=None, angle=None,
start_point=None, end_point=None, start_point=None, end_point=None,
latlon=False, cache=None, meta=True): latlon=False, cache=None, meta=True):
@ -165,6 +166,16 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
coordinates, and must be specified if *wrfin* is None. Default coordinates, and must be specified if *wrfin* is None. Default
is None. is None.
ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower
left latitude(s) for your domain, and must be specified if
*wrfin* is None. If the domain is a moving nest, this should be a
sequence of lower left latitudes. Default is None.
ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower
left longitude(s) for your domain, and must be specified if
*wrfin* is None. If the domain is a moving nest, this should be
a sequence of lower left longitudes. Default is None.
pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for
the pivot point, which indicates the location through which the pivot point, which indicates the location through which
the plane will pass. Must also specify *angle*. The coordinate the plane will pass. Must also specify *angle*. The coordinate
@ -291,6 +302,7 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
@set_interp_metadata("line") @set_interp_metadata("line")
def interpline(field2d, pivot_point=None, def interpline(field2d, pivot_point=None,
wrfin=None, timeidx=0, stagger=None, projection=None, wrfin=None, timeidx=0, stagger=None, projection=None,
ll_lat=None, ll_lon=None,
angle=None, start_point=None, angle=None, start_point=None,
end_point=None, latlon=False, end_point=None, latlon=False,
cache=None, meta=True): cache=None, meta=True):
@ -340,6 +352,16 @@ def interpline(field2d, pivot_point=None,
not be used when working with x,y coordinates. Default not be used when working with x,y coordinates. Default
is None. is None.
ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower
left latitude(s) for your domain, and must be specified if
*wrfin* is None. If the domain is a moving nest, this should be a
sequence of lower left latitudes. Default is None.
ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower
left longitude(s) for your domain, and must be specified if
*wrfin* is None. If the domain is a moving nest, this should be
a sequence of lower left longitudes. Default is None.
pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for
the pivot point, which indicates the location through which the pivot point, which indicates the location through which
the plane will pass. Must also specify *angle*. The coordinate the plane will pass. Must also specify *angle*. The coordinate
@ -420,7 +442,7 @@ def interpline(field2d, pivot_point=None,
if pivot_point is not None: if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None: if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
pivot_point_xy = (xy_coords.x, xy_coords.y) pivot_point_xy = (xy_coords.x, xy_coords.y)
else: else:
pivot_point_xy = (pivot_point.x, pivot_point.y) pivot_point_xy = (pivot_point.x, pivot_point.y)
@ -428,14 +450,14 @@ def interpline(field2d, pivot_point=None,
if start_point is not None and end_point is not None: if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None: if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx, xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
start_point_xy = (xy_coords.x, xy_coords.y) start_point_xy = (xy_coords.x, xy_coords.y)
else: else:
start_point_xy = (start_point.x, start_point.y) start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None: if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx, xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
end_point_xy = (xy_coords.x, xy_coords.y) end_point_xy = (xy_coords.x, xy_coords.y)
else: else:
end_point_xy = (end_point.x, end_point.y) end_point_xy = (end_point.x, end_point.y)

24
src/wrf/interputils.py
Unescape View File

@ -329,7 +329,8 @@ def get_xy(var, pivot_point=None, angle=None,
return xy return xy
def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None): def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None,
ll_lat=None, ll_lon=None):
"""Return the coordinate pairs in grid space. """Return the coordinate pairs in grid space.
This function converts latitude,longitude coordinate pairs to This function converts latitude,longitude coordinate pairs to
@ -371,6 +372,16 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None):
projection object to use when working with latitude, longitude projection object to use when working with latitude, longitude
coordinates, and must be specified if *wrfin* is None. Default coordinates, and must be specified if *wrfin* is None. Default
is None. is None.
ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower
left latitude(s) for your domain, and must be specified if
*wrfin* is None. If the domain is a moving nest, this should be a
sequence of lower left latitudes. Default is None.
ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower
left longitude(s) for your domain, and must be specified if
*wrfin* is None. If the domain is a moving nest, this should be
a sequence of lower left longitudes. Default is None.
Returns: Returns:
@ -379,8 +390,11 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None):
""" """
if wrfin is None and projection is None: if (wrfin is None and
raise ValueError ("'wrfin' or 'projection' parameter is required") (projection is None or ll_lat is None or ll_lon is None)):
raise ValueError ("'wrfin' parameter or "
"'projection', 'll_lat', and 'll_lon' parameters "
"are required")
if isinstance(pairs, Iterable): if isinstance(pairs, Iterable):
lat = [pair.lat for pair in pairs] lat = [pair.lat for pair in pairs]
@ -415,8 +429,8 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None):
truelat1=projection.truelat1, truelat1=projection.truelat1,
truelat2=projection.truelat2, truelat2=projection.truelat2,
stand_lon=projection.stand_lon, stand_lon=projection.stand_lon,
ref_lat=projection.ll_lat, ref_lat=ll_lat,
ref_lon=projection.ll_lon, ref_lon=ll_lon,
pole_lat=pole_lat, pole_lat=pole_lat,
pole_lon=pole_lon, pole_lon=pole_lon,
known_x=0, known_x=0,

18
src/wrf/metadecorators.py
Unescape View File

@ -887,6 +887,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
argvars = from_args(wrapped, ("field3d", "vert", "levels", argvars = from_args(wrapped, ("field3d", "vert", "levels",
"latlon", "missing", "latlon", "missing",
"wrfin", "timeidx", "stagger", "projection", "wrfin", "timeidx", "stagger", "projection",
"ll_lat", "ll_lon",
"pivot_point", "angle", "pivot_point", "angle",
"start_point", "end_point", "start_point", "end_point",
"cache"), "cache"),
@ -901,6 +902,8 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
timeidx = argvars["timeidx"] timeidx = argvars["timeidx"]
stagger = argvars["stagger"] stagger = argvars["stagger"]
projection = argvars["projection"] projection = argvars["projection"]
ll_lat = argvars["ll_lat"]
ll_lon = argvars["ll_lon"]
pivot_point = argvars["pivot_point"] pivot_point = argvars["pivot_point"]
angle = argvars["angle"] angle = argvars["angle"]
start_point = argvars["start_point"] start_point = argvars["start_point"]
@ -914,7 +917,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
if pivot_point is not None: if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None: if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
pivot_point_xy = (xy_coords.x, xy_coords.y) pivot_point_xy = (xy_coords.x, xy_coords.y)
else: else:
pivot_point_xy = (pivot_point.x, pivot_point.y) pivot_point_xy = (pivot_point.x, pivot_point.y)
@ -922,14 +925,14 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
if start_point is not None and end_point is not None: if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None: if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx, xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
start_point_xy = (xy_coords.x, xy_coords.y) start_point_xy = (xy_coords.x, xy_coords.y)
else: else:
start_point_xy = (start_point.x, start_point.y) start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None: if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx, xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
end_point_xy = (xy_coords.x, xy_coords.y) end_point_xy = (xy_coords.x, xy_coords.y)
else: else:
end_point_xy = (end_point.x, end_point.y) end_point_xy = (end_point.x, end_point.y)
@ -1104,6 +1107,7 @@ def _set_line_meta(wrapped, instance, args, kwargs):
""" """
argvars = from_args(wrapped, ("field2d", argvars = from_args(wrapped, ("field2d",
"wrfin", "timeidx", "stagger", "projection", "wrfin", "timeidx", "stagger", "projection",
"ll_lat", "ll_lon",
"pivot_point", "angle", "pivot_point", "angle",
"start_point", "end_point", "latlon", "start_point", "end_point", "latlon",
"cache"), "cache"),
@ -1114,6 +1118,8 @@ def _set_line_meta(wrapped, instance, args, kwargs):
timeidx = argvars["timeidx"] timeidx = argvars["timeidx"]
stagger = argvars["stagger"] stagger = argvars["stagger"]
projection = argvars["projection"] projection = argvars["projection"]
ll_lat = argvars["ll_lat"]
ll_lon = argvars["ll_lon"]
pivot_point = argvars["pivot_point"] pivot_point = argvars["pivot_point"]
angle = argvars["angle"] angle = argvars["angle"]
start_point = argvars["start_point"] start_point = argvars["start_point"]
@ -1131,7 +1137,7 @@ def _set_line_meta(wrapped, instance, args, kwargs):
if pivot_point is not None: if pivot_point is not None:
if pivot_point.lat is not None and pivot_point.lon is not None: if pivot_point.lat is not None and pivot_point.lon is not None:
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
pivot_point_xy = (xy_coords.x, xy_coords.y) pivot_point_xy = (xy_coords.x, xy_coords.y)
else: else:
pivot_point_xy = (pivot_point.x, pivot_point.y) pivot_point_xy = (pivot_point.x, pivot_point.y)
@ -1140,14 +1146,14 @@ def _set_line_meta(wrapped, instance, args, kwargs):
if start_point is not None and end_point is not None: if start_point is not None and end_point is not None:
if start_point.lat is not None and start_point.lon is not None: if start_point.lat is not None and start_point.lon is not None:
xy_coords = to_xy_coords(start_point, wrfin, timeidx, xy_coords = to_xy_coords(start_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
start_point_xy = (xy_coords.x, xy_coords.y) start_point_xy = (xy_coords.x, xy_coords.y)
else: else:
start_point_xy = (start_point.x, start_point.y) start_point_xy = (start_point.x, start_point.y)
if end_point.lat is not None and end_point.lon is not None: if end_point.lat is not None and end_point.lon is not None:
xy_coords = to_xy_coords(end_point, wrfin, timeidx, xy_coords = to_xy_coords(end_point, wrfin, timeidx,
stagger, projection) stagger, projection, ll_lat, ll_lon)
end_point_xy = (xy_coords.x, xy_coords.y) end_point_xy = (xy_coords.x, xy_coords.y)
else: else:
end_point_xy = (end_point.x, end_point.y) end_point_xy = (end_point.x, end_point.y)

32
test/utests.py
Unescape View File

@ -293,7 +293,23 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
# Note: Until the bug is fixed in NCL, the wrf-python cross # Note: Until the bug is fixed in NCL, the wrf-python cross
# sections will contain one extra point # sections will contain one extra point
nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross, rtol=.01) nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross,
rtol=.01)
# Test the manual projection method with lat/lon
lats = hts.coords["XLAT"]
lons = hts.coords["XLONG"]
ll_lat = lats[0,0]
ll_lon = lons[0,0]
pivot = CoordPair(lat=lats[lats.shape[-2]/2, lats.shape[-1]/2],
lon=lons[lons.shape[-2]/2, lons.shape[-1]/2])
v1 = vertcross(hts,p,wrfin=in_wrfnc,pivot_point=pivot_point,
angle=90.0)
v2 = vertcross(hts,p,projection=hts.attrs["projection"],
ll_lat=ll_lat, ll_lon=ll_lon,
pivot_point=pivot_point, angle=90.)
nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01)
# Test opposite # Test opposite
p_cross1 = vertcross(p,hts,pivot_point=pivot_point, angle=90.0) p_cross1 = vertcross(p,hts,pivot_point=pivot_point, angle=90.0)
@ -329,6 +345,20 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
# Note: After NCL is fixed, remove the slice. # Note: After NCL is fixed, remove the slice.
nt.assert_allclose(to_np(t2_line1)[...,0:-1], ref_t2_line) nt.assert_allclose(to_np(t2_line1)[...,0:-1], ref_t2_line)
# Test the manual projection method with lat/lon
lats = t2.coords["XLAT"]
lons = t2.coords["XLONG"]
ll_lat = lats[0,0]
ll_lon = lons[0,0]
pivot = CoordPair(lat=lats[lats.shape[-2]/2, lats.shape[-1]/2],
lon=lons[lons.shape[-2]/2, lons.shape[-1]/2])
l1 = interpline(t2,wrfin=in_wrfnc,pivot_point=pivot_point,
angle=90.0)
l2 = interpline(t2,projection=t2.attrs["projection"],
ll_lat=ll_lat, ll_lon=ll_lon,
pivot_point=pivot_point, angle=90.)
nt.assert_allclose(to_np(l1), to_np(l2), rtol=.01)
# Test point to point # Test point to point
start_point = CoordPair(0, t2.shape[-2]/2) start_point = CoordPair(0, t2.shape[-2]/2)
end_point = CoordPair(-1, t2.shape[-2]/2) end_point = CoordPair(-1, t2.shape[-2]/2)

Write Preview
Loading…
Cancel
Save