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

Added autolevels parameter for vertcross 

Users can now specify how many autogenerated levels they want

Updated unit tests with NCL for the new routines.
lon0
Bill Ladwig 7 years ago
parent
70371d191b
commit
a3880b1e10
5 changed files with 146 additions and 44 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. 9
      src/wrf/interp.py
  2. 18
      src/wrf/interputils.py
  3. 5
      src/wrf/metadecorators.py
  4. 37
      test/ncl_get_var.ncl
  5. 121
      test/utests.py

9
src/wrf/interp.py
Unescape View File

@ -94,7 +94,7 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
ll_point=None, ll_point=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, autolevels=100, cache=None, meta=True):
"""Return the vertical cross section for a three-dimensional field. """Return the vertical cross section for a three-dimensional field.
The cross section is defined by a horizontal line through the domain. The cross section is defined by a horizontal line through the domain.
@ -214,6 +214,10 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
and the latitude and longitude information will not be and the latitude and longitude information will not be
present. present.
autolevels(:obj:`int`, optional): The number of evenly spaced
automatically chosen vertical levels to use when *levels*
is None. Default is 100.
cache (:obj:`dict`, optional): A dictionary of (varname, ndarray) cache (:obj:`dict`, optional): A dictionary of (varname, ndarray)
that can be used to supply pre-extracted NetCDF variables to the that can be used to supply pre-extracted NetCDF variables to the
computational routines. It is primarily used for internal computational routines. It is primarily used for internal
@ -279,7 +283,8 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
xy, var2dz, z_var2d = get_xy_z_params(to_np(vert), pivot_point_xy, xy, var2dz, z_var2d = get_xy_z_params(to_np(vert), pivot_point_xy,
angle, start_point_xy, angle, start_point_xy,
end_point_xy, levels) end_point_xy, levels,
autolevels)
if not multi: if not multi:
result = _vertcross(field3d, xy, var2dz, z_var2d, missing) result = _vertcross(field3d, xy, var2dz, z_var2d, missing)

18
src/wrf/interputils.py
Unescape View File

@ -185,7 +185,7 @@ def _calc_xy(xdim, ydim, pivot_point=None, angle=None,
def get_xy_z_params(z, pivot_point=None, angle=None, def get_xy_z_params(z, pivot_point=None, angle=None,
start_point=None, end_point=None, start_point=None, end_point=None,
levels=None): levels=None, autolevels=100):
"""Return the cross section parameters. """Return the cross section parameters.
This function returns the xy horizontal cross section line coordinates, This function returns the xy horizontal cross section line coordinates,
@ -224,6 +224,10 @@ def get_xy_z_params(z, pivot_point=None, angle=None,
vertical levels in the output array. If None, a fixed set of vertical levels in the output array. If None, a fixed set of
vertical levels is provided. Default is None. vertical levels is provided. Default is None.
autolevels(:obj:`int`, optional): The number of evenly spaced
automatically chosen vertical levels to use when *levels*
is None. Default is 100.
Returns: Returns:
:obj:`tuple`: A tuple containing the xy horizontal cross section :obj:`tuple`: A tuple containing the xy horizontal cross section
@ -248,20 +252,18 @@ def get_xy_z_params(z, pivot_point=None, angle=None,
if(var2dz[idx1] > var2dz[idx2]): # monotonically decreasing coordinate if(var2dz[idx1] > var2dz[idx2]): # monotonically decreasing coordinate
z_max = floor(np.amax(z)/10) * 10 # bottom value z_max = floor(np.amax(z)/10) * 10 # bottom value
z_min = ceil(np.amin(z)/10) * 10 # top value z_min = ceil(np.amin(z)/10) * 10 # top value
dz = 10 dz = (1.0/autolevels) * (z_max - z_min)
nlevels = int((z_max - z_min)/dz) z_var2d = np.zeros((autolevels), dtype=z.dtype)
z_var2d = np.zeros((nlevels), dtype=z.dtype)
z_var2d[0] = z_max z_var2d[0] = z_max
dz = -dz dz = -dz
else: else:
z_max = np.amax(z) z_max = np.amax(z)
z_min = 0. z_min = 0.
dz = 0.01*z_max dz = (1.0/autolevels)*z_max
nlevels = int(z_max/dz) z_var2d = np.zeros((autolevels), dtype=z.dtype)
z_var2d = np.zeros((nlevels), dtype=z.dtype)
z_var2d[0] = z_min z_var2d[0] = z_min
for i in py3range(1,nlevels): for i in py3range(1,autolevels):
z_var2d[i] = z_var2d[0] + i*dz z_var2d[i] = z_var2d[0] + i*dz
else: else:
z_var2d = np.asarray(levels, z.dtype) z_var2d = np.asarray(levels, z.dtype)

5
src/wrf/metadecorators.py
Unescape View File

@ -889,7 +889,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
"latlon", "missing", "latlon", "missing",
"wrfin", "timeidx", "stagger", "projection", "wrfin", "timeidx", "stagger", "projection",
"ll_point", "pivot_point", "angle", "ll_point", "pivot_point", "angle",
"start_point", "end_point", "start_point", "end_point", "autolevels",
"cache"), "cache"),
*args, **kwargs) *args, **kwargs)
@ -907,6 +907,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
angle = argvars["angle"] angle = argvars["angle"]
start_point = argvars["start_point"] start_point = argvars["start_point"]
end_point = argvars["end_point"] end_point = argvars["end_point"]
autolevels = argvars["autolevels"]
cache = argvars["cache"] cache = argvars["cache"]
start_point_xy = None start_point_xy = None
@ -937,7 +938,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
end_point_xy = (end_point.x, end_point.y) end_point_xy = (end_point.x, end_point.y)
xy, var2dz, z_var2d = get_xy_z_params(to_np(z), pivot_point_xy, angle, xy, var2dz, z_var2d = get_xy_z_params(to_np(z), pivot_point_xy, angle,
start_point_xy, end_point_xy, levels) start_point_xy, end_point_xy, levels, autolevels)
# Make a copy so we don't modify a user supplied cache # Make a copy so we don't modify a user supplied cache
if cache is not None: if cache is not None:

37
test/ncl_get_var.ncl
Unescape View File

@ -5,7 +5,7 @@
;system("printenv") ;system("printenv")
if (.not. isvar("in_file")) then if (.not. isvar("in_file")) then
in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d02_2010-06-13_21:00:00.nc" in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc"
end if end if
if (.not. isvar("out_file")) then if (.not. isvar("out_file")) then
@ -48,20 +48,21 @@
; Do the interpolation routines manually ; Do the interpolation routines manually
; 3D vertical cross section ; 3D vertical cross section
z = wrf_user_getvar(input_file, "z", 0) ; grid point height z := wrf_user_getvar(input_file, "z", 0) ; grid point height
p = wrf_user_getvar(input_file, "pressure", 0) ; total pressure p := wrf_user_getvar(input_file, "pressure", 0) ; total pressure
dimsz = dimsizes(z) dimsz = dimsizes(z)
pivot = (/ dimsz(2)/2, dimsz(1)/2 /) ; pivot point is center of domain pivot = (/ dimsz(2)/2, dimsz(1)/2 /) ; pivot point is center of domain
ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False) ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False)
p_cross = wrf_user_intrp3d(p,z,"v",pivot,90.0,False) p_cross = wrf_user_intrp3d(p,z,"v",pivot,90.0,False)
p_cross!0 = "Vertical_p" p_cross!0 = "Vertical_p"
p_cross!1 = "Horizontal_p" p_cross!1 = "Horizontal_p"
fout->ht_cross = ht_cross fout->ht_cross = ht_cross
fout->p_cross = p_cross fout->p_cross = p_cross
;
; For the new cross section routine ; For the new cross section routine
xopt = True xopt = True
xopt@use_pivot = True xopt@use_pivot = True
@ -90,17 +91,17 @@
fout->ht_vertcross2 = ht_vertcross2 fout->ht_vertcross2 = ht_vertcross2
; Note for the default file: lats = wrf_user_getvar(input_file, "lat", 0)
; >>> np.amin(lats) lons = wrf_user_getvar(input_file, "lon", 0)
; 32.77895
; >>> np.amax(lats) start_lat = min(lats) + .25d*(max(lats) - min(lats))
; 40.209736 end_lat = min(lats) + .75d*(max(lats) - min(lats))
; >>> np.amin(lons)
; -104.71738 start_lon = min(lons) + .25d*(max(lons) - min(lons))
; >>> np.amax(lons) end_lon = min(lons) + .75d*(max(lons) - min(lons))
; -95.153076
start_end = (/ start_lon, start_lat, end_lon, end_lat /)
start_end = (/ -102.0, 35.0, -98.0, 37.0 /)
; For the new cross section routine ; For the new cross section routine
xopt := True xopt := True
xopt@use_pivot = False xopt@use_pivot = False
@ -108,6 +109,7 @@
xopt@file_handle = input_file xopt@file_handle = input_file
xopt@timeidx = 0 xopt@timeidx = 0
xopt@linecoords = True xopt@linecoords = True
xopt@autolevels = 1000
ht_vertcross3 = wrf_user_vertcross(z, p, start_end, xopt) ht_vertcross3 = wrf_user_vertcross(z, p, start_end, xopt)
ht_vertcross3!0 = "vertical3" ht_vertcross3!0 = "vertical3"
@ -118,7 +120,6 @@
; 3D horizontal interpolation ; 3D horizontal interpolation
plev = 500. ; 500 MB plev = 500. ; 500 MB
z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False) z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False)
;z_500_dims = dimsizes(z_500)
fout->z_500 = z_500 fout->z_500 = z_500
@ -221,8 +222,8 @@
; lat/lon to x/y and x/y to lat/lon routines ; lat/lon to x/y and x/y to lat/lon routines
lats = (/-55, -60, -65 /) lats := (/-55, -60, -65 /)
lons = (/25, 30, 35 /) lons := (/25, 30, 35 /)
i_s = (/10, 100, 150 /) i_s = (/10, 100, 150 /)
j_s = (/10, 100, 150 /) j_s = (/10, 100, 150 /)
ij = wrf_user_ll_to_ij(input_file, lons, lats, True) ij = wrf_user_ll_to_ij(input_file, lons, lats, True)

121
test/utests.py
Unescape View File

@ -65,6 +65,10 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
timeidx = 0 timeidx = 0
if not pynio: if not pynio:
in_wrfnc = NetCDF(wrf_in) in_wrfnc = NetCDF(wrf_in)
try:
in_wrfnc.set_auto_mask(False)
except:
pass
else: else:
# Note: Python doesn't like it if you reassign an outer scoped # Note: Python doesn't like it if you reassign an outer scoped
# variable (wrf_in in this case) # variable (wrf_in in this case)
@ -77,6 +81,10 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
timeidx = None timeidx = None
if not pynio: if not pynio:
nc = NetCDF(wrf_in) nc = NetCDF(wrf_in)
try:
nc.set_auto_mask(False)
except:
pass
in_wrfnc = [nc for i in xrange(repeat)] in_wrfnc = [nc for i in xrange(repeat)]
else: else:
if not wrf_in.endswith(".nc"): if not wrf_in.endswith(".nc"):
@ -88,6 +96,10 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
refnc = NetCDF(referent) refnc = NetCDF(referent)
try:
refnc.set_auto_mask(False)
except:
pass
# These have a left index that defines the product type # These have a left index that defines the product type
multiproduct = varname in ("uvmet", "uvmet10", "cape_2d", "cape_3d", multiproduct = varname in ("uvmet", "uvmet10", "cape_2d", "cape_3d",
@ -204,6 +216,10 @@ def _get_refvals(referent, varname, repeat, multi):
pass pass
refnc = NetCDF(referent) refnc = NetCDF(referent)
try:
refnc.set_auto_mask(False)
except:
pass
if not multi: if not multi:
ref_vals = refnc.variables[varname][:] ref_vals = refnc.variables[varname][:]
@ -244,6 +260,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
timeidx = 0 timeidx = 0
if not pynio: if not pynio:
in_wrfnc = NetCDF(wrf_in) in_wrfnc = NetCDF(wrf_in)
try:
in_wrfnc.set_auto_mask(False)
except:
pass
else: else:
# Note: Python doesn't like it if you reassign an outer scoped # Note: Python doesn't like it if you reassign an outer scoped
# variable (wrf_in in this case) # variable (wrf_in in this case)
@ -256,6 +276,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
timeidx = None timeidx = None
if not pynio: if not pynio:
nc = NetCDF(wrf_in) nc = NetCDF(wrf_in)
try:
nc.set_auto_mask(False)
except:
pass
in_wrfnc = [nc for i in xrange(repeat)] in_wrfnc = [nc for i in xrange(repeat)]
else: else:
if not wrf_in.endswith(".nc"): if not wrf_in.endswith(".nc"):
@ -279,22 +303,35 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
elif (varname == "vertcross"): elif (varname == "vertcross"):
ref_ht_cross = _get_refvals(referent, "ht_cross", repeat, multi) ref_ht_cross = _get_refvals(referent, "ht_cross", repeat, multi)
ref_p_cross = _get_refvals(referent, "p_cross", repeat, multi) ref_p_cross = _get_refvals(referent, "p_cross", repeat, multi)
ref_ht_vertcross1 = _get_refvals(referent, "ht_vertcross1", repeat,
multi)
ref_ht_vertcross2 = _get_refvals(referent, "ht_vertcross2", repeat,
multi)
ref_ht_vertcross3 = _get_refvals(referent, "ht_vertcross3", repeat,
multi)
hts = getvar(in_wrfnc, "z", timeidx=timeidx) hts = getvar(in_wrfnc, "z", timeidx=timeidx)
p = getvar(in_wrfnc, "pressure", timeidx=timeidx) p = getvar(in_wrfnc, "pressure", timeidx=timeidx)
pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2) pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2)
#ht_cross = vertcross(to_np(hts), p, pivot_point=pivot_point,
# angle=90., latlon=True) # Beginning in wrf-python 1.3, users can select number of levels.
# Originally, for pressure, dz was 10, so let's back calculate
# the number of levels.
p_max = np.floor(np.amax(p)/10) * 10 # bottom value
p_min = np.ceil(np.amin(p)/10) * 10 # top value
p_autolevels = int((p_max - p_min) /10)
# Make sure the numpy versions work first # Make sure the numpy versions work first
ht_cross = vertcross(to_np(hts), to_np(p), ht_cross = vertcross(to_np(hts), to_np(p),
pivot_point=pivot_point, angle=90.) pivot_point=pivot_point, angle=90.,
ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.) autolevels=p_autolevels)
ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.,
autolevels=p_autolevels)
# Note: Until the bug is fixed in NCL, the wrf-python cross nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01)
# sections will contain one extra point
nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross,
rtol=.01)
# Test the manual projection method with lat/lon # Test the manual projection method with lat/lon
lats = hts.coords["XLAT"] lats = hts.coords["XLAT"]
@ -304,6 +341,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
int(lats.shape[-1]/2)], int(lats.shape[-1]/2)],
lon=lons[int(lons.shape[-2]/2), lon=lons[int(lons.shape[-2]/2),
int(lons.shape[-1]/2)]) int(lons.shape[-1]/2)])
v1 = vertcross(hts,p,wrfin=in_wrfnc,pivot_point=pivot_point, v1 = vertcross(hts,p,wrfin=in_wrfnc,pivot_point=pivot_point,
angle=90.0) angle=90.0)
v2 = vertcross(hts,p,projection=hts.attrs["projection"], v2 = vertcross(hts,p,projection=hts.attrs["projection"],
@ -313,9 +351,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01) 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)
nt.assert_allclose(to_np(p_cross1)[...,0:-1], nt.assert_allclose(to_np(p_cross1),
ref_p_cross, ref_p_cross,
rtol=.01) rtol=.01)
# Test point to point # Test point to point
@ -329,6 +368,50 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
nt.assert_allclose(to_np(p_cross1), nt.assert_allclose(to_np(p_cross1),
to_np(p_cross2)) to_np(p_cross2))
# Check the new vertcross routine
pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2)
ht_cross = vertcross(hts, p,
pivot_point=pivot_point, angle=90.,
latlon=True)
nt.assert_allclose(to_np(ht_cross),
to_np(ref_ht_vertcross1), atol=.01)
levels = [1000., 850., 700., 500., 250.]
ht_cross = vertcross(hts, p,
pivot_point=pivot_point, angle=90.,
levels=levels, latlon=True)
nt.assert_allclose(to_np(ht_cross),
to_np(ref_ht_vertcross2), atol=.01)
start_lat = np.amin(lats) + .25*(np.amax(lats) - np.amin(lats))
end_lat = np.amin(lats) + .75*(np.amax(lats) - np.amin(lats))
start_lon = np.amin(lons) + .25*(np.amax(lons) - np.amin(lons))
end_lon = np.amin(lons) + .75*(np.amax(lons) - np.amin(lons))
start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)
# ll_point and projection came from above
ht_cross = vertcross(hts, p,
start_point=start_point,
end_point=end_point,
projection=hts.attrs["projection"],
ll_point=ll_point,
latlon=True,
autolevels=1000)
nt.assert_allclose(to_np(ht_cross),
to_np(ref_ht_vertcross3),
rtol=.01)
elif (varname == "interpline"): elif (varname == "interpline"):
ref_t2_line = _get_refvals(referent, "t2_line", repeat, multi) ref_t2_line = _get_refvals(referent, "t2_line", repeat, multi)
@ -336,15 +419,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
t2 = getvar(in_wrfnc, "T2", timeidx=timeidx) t2 = getvar(in_wrfnc, "T2", timeidx=timeidx)
pivot_point = CoordPair(t2.shape[-1] / 2, t2.shape[-2] / 2) pivot_point = CoordPair(t2.shape[-1] / 2, t2.shape[-2] / 2)
#t2_line1 = interpline(to_np(t2), pivot_point=pivot_point,
# angle=90.0, latlon=True)
# Make sure the numpy version works # Make sure the numpy version works
t2_line1 = interpline(to_np(t2), pivot_point=pivot_point, t2_line1 = interpline(to_np(t2), pivot_point=pivot_point,
angle=90.0) angle=90.0)
t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0) t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0)
# Note: After NCL is fixed, remove the slice. nt.assert_allclose(to_np(t2_line1), ref_t2_line)
nt.assert_allclose(to_np(t2_line1)[...,0:-1], ref_t2_line)
# Test the manual projection method with lat/lon # Test the manual projection method with lat/lon
lats = t2.coords["XLAT"] lats = t2.coords["XLAT"]
@ -575,6 +655,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
timeidx = 0 timeidx = 0
if not pynio: if not pynio:
in_wrfnc = NetCDF(wrf_in) in_wrfnc = NetCDF(wrf_in)
try:
in_wrfnc.set_auto_mask(False)
except:
pass
else: else:
# Note: Python doesn't like it if you reassign an outer scoped # Note: Python doesn't like it if you reassign an outer scoped
# variable (wrf_in in this case) # variable (wrf_in in this case)
@ -587,6 +671,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
timeidx = None timeidx = None
if not pynio: if not pynio:
nc = NetCDF(wrf_in) nc = NetCDF(wrf_in)
try:
nc.set_auto_mask(False)
except:
pass
in_wrfnc = [nc for i in xrange(repeat)] in_wrfnc = [nc for i in xrange(repeat)]
else: else:
if not wrf_in.endswith(".nc"): if not wrf_in.endswith(".nc"):
@ -597,6 +685,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
in_wrfnc = [nc for i in xrange(repeat)] in_wrfnc = [nc for i in xrange(repeat)]
refnc = NetCDF(referent) refnc = NetCDF(referent)
try:
refnc.set_auto_mask(False)
except:
pass
if testid == "xy": if testid == "xy":
# Since this domain is not moving, the reference values are the # Since this domain is not moving, the reference values are the
@ -691,7 +783,8 @@ if __name__ == "__main__":
"geopt", "helicity", "lat", "lon", "omg", "p", "pressure", "geopt", "helicity", "lat", "lon", "omg", "p", "pressure",
"pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc",
"theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va",
"wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"] "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag",
"height_agl"]
interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"] interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"]
latlon_tests = ["xy", "ll"] latlon_tests = ["xy", "ll"]

Write Preview
Loading…
Cancel
Save