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

Fixed metadata bugs related to interpolation. 

Unit tests were also updated to test for this situation.
lon0
Bill Ladwig 7 years ago
parent
5d448aa278
commit
0d6f9230da
5 changed files with 58 additions and 26 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. 32
      src/wrf/metadecorators.py
  2. BIN
      test/ci_tests/ci_test_file.nc
  3. 6
      test/ci_tests/make_test_file.py
  4. 23
      test/ci_tests/utests.py
  5. 23
      test/utests.py

32
src/wrf/metadecorators.py
Unescape View File

@ -805,7 +805,7 @@ def _set_horiz_meta(wrapped, instance, args, kwargs):
outname = None outname = None
outdimnames = None outdimnames = None
outcoords = None outcoords = None
outattrs = None outattrs = OrderedDict()
# Get the vertical level units # Get the vertical level units
vert_units = None vert_units = None
@ -823,7 +823,6 @@ def _set_horiz_meta(wrapped, instance, args, kwargs):
if isinstance(field3d, DataArray): if isinstance(field3d, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outattrs = OrderedDict()
outdimnames = list(field3d.dims) outdimnames = list(field3d.dims)
outcoords.update(field3d.coords) outcoords.update(field3d.coords)
outdimnames.remove(field3d.dims[-3]) outdimnames.remove(field3d.dims[-3])
@ -836,7 +835,6 @@ def _set_horiz_meta(wrapped, instance, args, kwargs):
else: else:
outname = "field3d_{0}".format(name_levelstr) outname = "field3d_{0}".format(name_levelstr)
outattrs = OrderedDict()
outattrs["level"] = levelstr outattrs["level"] = levelstr
outattrs["missing_value"] = missingval outattrs["missing_value"] = missingval
@ -959,7 +957,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
outname = None outname = None
outdimnames = None outdimnames = None
outcoords = None outcoords = None
outattrs = None outattrs = OrderedDict()
# Use XY to set the cross-section metadata # Use XY to set the cross-section metadata
st_x = xy[0,0] st_x = xy[0,0]
@ -974,7 +972,6 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
if isinstance(field3d, DataArray): if isinstance(field3d, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outattrs = OrderedDict()
outdimnames = list(field3d.dims) outdimnames = list(field3d.dims)
outcoords.update(field3d.coords) outcoords.update(field3d.coords)
for i in py3range(-3,0,1): for i in py3range(-3,0,1):
@ -1062,7 +1059,6 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
"not of type xarray.DataArray and contains no " "not of type xarray.DataArray and contains no "
"coordinate information") "coordinate information")
outname = "field3d_cross" outname = "field3d_cross"
outattrs = OrderedDict()
outattrs["orientation"] = cross_str outattrs["orientation"] = cross_str
outattrs["missing_value"] = missingval outattrs["missing_value"] = missingval
@ -1173,7 +1169,7 @@ def _set_line_meta(wrapped, instance, args, kwargs):
outname = None outname = None
outdimnames = None outdimnames = None
outcoords = None outcoords = None
outattrs = None outattrs = OrderedDict()
# Use XY to set the cross-section metadata # Use XY to set the cross-section metadata
st_x = xy[0,0] st_x = xy[0,0]
@ -1188,7 +1184,6 @@ def _set_line_meta(wrapped, instance, args, kwargs):
if isinstance(field2d, DataArray): if isinstance(field2d, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outattrs = OrderedDict()
outdimnames = list(field2d.dims) outdimnames = list(field2d.dims)
outcoords.update(field2d.coords) outcoords.update(field2d.coords)
for i in py3range(-2,0,1): for i in py3range(-2,0,1):
@ -1273,7 +1268,6 @@ def _set_line_meta(wrapped, instance, args, kwargs):
"not of type xarray.DataArray and contains no " "not of type xarray.DataArray and contains no "
"coordinate information") "coordinate information")
outname = "field2d_line" outname = "field2d_line"
outattrs = OrderedDict()
outattrs["orientation"] = cross_str outattrs["orientation"] = cross_str
@ -1334,12 +1328,11 @@ def _set_vinterp_meta(wrapped, instance, args, kwargs):
outname = None outname = None
outdimnames = None outdimnames = None
outcoords = None outcoords = None
outattrs = None outattrs = OrderedDict()
if isinstance(field, DataArray): if isinstance(field, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outattrs = OrderedDict()
outdimnames = list(field.dims) outdimnames = list(field.dims)
outcoords.update(field.coords) outcoords.update(field.coords)
@ -1352,13 +1345,14 @@ def _set_vinterp_meta(wrapped, instance, args, kwargs):
outdimnames.insert(-2, "interp_level") outdimnames.insert(-2, "interp_level")
outcoords["interp_level"] = interp_levels outcoords["interp_level"] = interp_levels
outattrs.update(field.attrs) outattrs.update(field.attrs)
outattrs["vert_interp_type"] = vert_coord
outname = field.name outname = field.name
else: else:
outname = field_type outname = field_type
outattrs["vert_interp_type"] = vert_coord
return DataArray(result, name=outname, dims=outdimnames, return DataArray(result, name=outname, dims=outdimnames,
coords=outcoords, attrs=outattrs) coords=outcoords, attrs=outattrs)
@ -1402,8 +1396,7 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs):
argvars = from_args(wrapped, ("field3d", "xy"), *args, **kwargs) argvars = from_args(wrapped, ("field3d", "xy"), *args, **kwargs)
field3d = argvars["field3d"] field3d = argvars["field3d"]
xy = argvars["xy"] xy = to_np(argvars["xy"])
xy = to_np(xy)
result = wrapped(*args, **kwargs) result = wrapped(*args, **kwargs)
@ -1419,12 +1412,11 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs):
outname = None outname = None
outdimnames = None outdimnames = None
outcoords = None outcoords = None
outattrs = None outattrs = OrderedDict()
# Dims are (...,xy,z) # Dims are (...,xy,z)
if isinstance(field3d, DataArray): if isinstance(field3d, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outattrs = OrderedDict()
outdimnames = list(field3d.dims) outdimnames = list(field3d.dims)
outcoords.update(field3d.coords) outcoords.update(field3d.coords)
@ -1520,12 +1512,11 @@ def _set_1d_meta(wrapped, instance, args, kwargs):
outname = None outname = None
outdimnames = None outdimnames = None
outcoords = None outcoords = None
outattrs = None outattrs = OrderedDict()
# Dims are (...,xy,z) # Dims are (...,xy,z)
if isinstance(field, DataArray): if isinstance(field, DataArray):
outcoords = OrderedDict() outcoords = OrderedDict()
outattrs = OrderedDict()
outdimnames = list(field.dims) outdimnames = list(field.dims)
outdimnames.pop(-1) outdimnames.pop(-1)
@ -1540,9 +1531,6 @@ def _set_1d_meta(wrapped, instance, args, kwargs):
outname = "{0}_z".format(field.name) outname = "{0}_z".format(field.name)
outcoords["z"] = z_out outcoords["z"] = z_out
outattrs["_FillValue"] = missingval
outattrs["missing_value"] = missingval
desc = field.attrs.get("description", None) desc = field.attrs.get("description", None)
if desc is not None: if desc is not None:
outattrs["description"] = desc outattrs["description"] = desc
@ -1554,6 +1542,8 @@ def _set_1d_meta(wrapped, instance, args, kwargs):
else: else:
outname = "field_z" outname = "field_z"
outattrs["_FillValue"] = missingval
outattrs["missing_value"] = missingval
return DataArray(result, name=outname, dims=outdimnames, return DataArray(result, name=outname, dims=outdimnames,
coords=outcoords, attrs=outattrs) coords=outcoords, attrs=outattrs)

BIN
test/ci_tests/ci_test_file.nc
View File

Binary file not shown.

6
test/ci_tests/make_test_file.py
Unescape View File

@ -9,9 +9,9 @@ from netCDF4 import Dataset
from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round, from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round,
CoordPair, ll_to_xy, xy_to_ll, to_np) CoordPair, ll_to_xy, xy_to_ll, to_np)
VARS_TO_KEEP = ("XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", "XLONG_V", VARS_TO_KEEP = ("Times", "XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U",
"U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", "T2", "XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2",
"PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD", "T2", "PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD",
"QGRAUP", "QRAIN", "QSNOW", "MAPFAC_M", "MAPFAC_U", "QGRAUP", "QRAIN", "QSNOW", "MAPFAC_M", "MAPFAC_U",
"MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC") "MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC")

23
test/ci_tests/utests.py
Unescape View File

@ -100,6 +100,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
ref_ht_850 = _get_refvals(referent, "interplevel", repeat, multi) ref_ht_850 = _get_refvals(referent, "interplevel", 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)
# Check that it works with numpy arrays
hts_850 = interplevel(to_np(hts), p, 850)
#print (hts_850)
hts_850 = interplevel(hts, p, 850) hts_850 = interplevel(hts, p, 850)
nt.assert_allclose(to_np(hts_850), ref_ht_850) nt.assert_allclose(to_np(hts_850), ref_ht_850)
@ -111,6 +115,10 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
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)
# Check that it works with numpy arrays
ht_cross = vertcross(to_np(hts), to_np(p),
pivot_point=pivot_point, angle=90.)
#print (ht_cross)
ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.) ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.)
nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01) nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01)
@ -123,6 +131,10 @@ 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)
# Check that it works with numpy arrays
t2_line1 = interpline(to_np(t2), pivot_point=pivot_point,
angle=90.0)
#print (t2_line1)
t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0) t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0)
nt.assert_allclose(to_np(t2_line1), ref_t2_line) nt.assert_allclose(to_np(t2_line1), ref_t2_line)
@ -136,6 +148,17 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
interp_levels = [200,300,500,1000] interp_levels = [200,300,500,1000]
# Check that it works with numpy arrays
field = vinterp(in_wrfnc,
field=to_np(tk),
vert_coord="theta",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
timeidx=timeidx,
log_p=True)
#print (field)
field = vinterp(in_wrfnc, field = vinterp(in_wrfnc,
field=tk, field=tk,
vert_coord="theta", vert_coord="theta",

23
test/utests.py
Unescape View File

@ -269,6 +269,9 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
ref_ht_500 = _get_refvals(referent, "z_500", repeat, multi) ref_ht_500 = _get_refvals(referent, "z_500", 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)
# Make sure the numpy versions work first
hts_500 = interplevel(to_np(hts), to_np(p), 500)
hts_500 = interplevel(hts, p, 500) hts_500 = interplevel(hts, p, 500)
nt.assert_allclose(to_np(hts_500), ref_ht_500) nt.assert_allclose(to_np(hts_500), ref_ht_500)
@ -283,6 +286,9 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
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, #ht_cross = vertcross(to_np(hts), p, pivot_point=pivot_point,
# angle=90., latlon=True) # angle=90., latlon=True)
# Make sure the numpy versions work first
ht_cross = vertcross(to_np(hts), to_np(p),
pivot_point=pivot_point, angle=90.)
ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.) ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.)
# Note: Until the bug is fixed in NCL, the wrf-python cross # Note: Until the bug is fixed in NCL, the wrf-python cross
@ -315,6 +321,9 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
#t2_line1 = interpline(to_np(t2), pivot_point=pivot_point, #t2_line1 = interpline(to_np(t2), pivot_point=pivot_point,
# angle=90.0, latlon=True) # angle=90.0, latlon=True)
# Make sure the numpy version works
t2_line1 = interpline(to_np(t2), pivot_point=pivot_point,
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. # Note: After NCL is fixed, remove the slice.
@ -337,6 +346,16 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
interp_levels = [200,300,500,1000] interp_levels = [200,300,500,1000]
# Make sure the numpy version works
field = vinterp(in_wrfnc,
field=to_np(tk),
vert_coord="theta",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
timeidx=timeidx,
log_p=True)
field = vinterp(in_wrfnc, field = vinterp(in_wrfnc,
field=tk, field=tk,
vert_coord="theta", vert_coord="theta",
@ -630,8 +649,8 @@ class WRFLatLonTest(ut.TestCase):
if __name__ == "__main__": if __name__ == "__main__":
from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule,
omp_set_dynamic, OMP_SCHED_STATIC) omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC)
omp_set_num_threads(8) omp_set_num_threads(omp_get_num_procs()-1)
omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_schedule(OMP_SCHED_STATIC, 0)
omp_set_dynamic(False) omp_set_dynamic(False)

Write Preview
Loading…
Cancel
Save