Skip to content

Commit 4eeeadc

Browse files
author
Bill Ladwig
committed
Fix issues with moving nests and line interpolation.
Updated unit tests.
1 parent c8ba485 commit 4eeeadc

4 files changed

Lines changed: 208 additions & 114 deletions

File tree

‎src/wrf/interp.py‎

Lines changed: 53 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -106,52 +106,6 @@ def interplevel(field3d, vert, desiredlev, missing=default_fill(np.float64),
106106
return masked
107107

108108

109-
def _vertcross_nest_alltimes(field3d, vert, levels=None,
110-
missing=default_fill(np.float64),
111-
wrfin=None, timeidx=0, stagger=None, projection=None,
112-
ll_point=None,
113-
pivot_point=None, angle=None,
114-
start_point=None, end_point=None,
115-
latlon=False, autolevels=100, cache=None, meta=True):
116-
117-
# Some fields like uvmet have an extra left dimension for the product
118-
# type, we'll handle that iteration here.
119-
multi = True if field3d.ndim - vert.ndim == 1 else False
120-
121-
# Check if we have a wrfin file, or this is a no go.
122-
if wrfin is None:
123-
raise ValueError("'wrfin' is required when using all times "
124-
"from a moving nest with lat/lon coords")
125-
126-
if multi:
127-
if field3d.ndim == 4:
128-
raise ValueError("all times requested for a moving nest, "
129-
"but no time dimension found for "
130-
'field3d')
131-
else:
132-
if field3d.ndim < 4:
133-
raise ValueError("all times requested for a moving nest, "
134-
"but no time dimension found for "
135-
'field3d')
136-
137-
numtimes = field3d.shape[-4]
138-
139-
for t in py3range(numtimes):
140-
#_meta = True if t == 0 else False
141-
_meta = True
142-
143-
v = vertcross(field3d, vert, levels, missing, wrfin, t, stagger,
144-
projection, ll_point, pivot_point, angle, start_point,
145-
end_point, latlon, autolevels, cache, _meta)
146-
147-
print(v.attrs)
148-
149-
150-
151-
152-
153-
154-
155109
@set_interp_metadata("cross")
156110
def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
157111
wrfin=None, timeidx=0, stagger=None, projection=None,
@@ -308,33 +262,6 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
308262
# type, we'll handle that iteration here.
309263
multi = True if field3d.ndim - vert.ndim == 1 else False
310264

311-
312-
if timeidx is None:
313-
if (latlon is True or is_latlon_pair(start_point) or
314-
is_latlon_pair(pivot_point)):
315-
316-
if wrfin is not None:
317-
# Moving nests aren't supported with ALL_TIMES because the
318-
# domain could move outside of the cross section, which causes
319-
# crashes or different line lengths.
320-
if is_moving_domain(wrfin):
321-
raise ValueError("Requesting all times with a moving nest "
322-
"is not supported when using lat/lon "
323-
"cross sections because the domain could "
324-
"move outside of the cross section. "
325-
"You must request each time "
326-
"individually.")
327-
else:
328-
_timeidx = 0
329-
330-
# If using grid coordinates, then don't care about lat/lon
331-
# coordinates. Just use 0.
332-
else:
333-
_timeidx = 0
334-
else:
335-
_timeidx = timeidx
336-
337-
338265
try:
339266
xy = cache["xy"]
340267
var2dz = cache["var2dz"]
@@ -345,6 +272,31 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
345272
end_point_xy = None
346273
pivot_point_xy = None
347274

275+
if timeidx is None:
276+
if (latlon is True or is_latlon_pair(start_point) or
277+
is_latlon_pair(pivot_point)):
278+
279+
if wrfin is not None:
280+
# Moving nests aren't supported with ALL_TIMES because the
281+
# domain could move outside of the cross section, which causes
282+
# crashes or different line lengths.
283+
if is_moving_domain(wrfin):
284+
raise ValueError("Requesting all times with a moving nest "
285+
"is not supported when using lat/lon "
286+
"cross sections because the domain could "
287+
"move outside of the cross section. "
288+
"You must request each time "
289+
"individually.")
290+
else:
291+
_timeidx = 0
292+
293+
# If using grid coordinates, then don't care about lat/lon
294+
# coordinates. Just use 0.
295+
else:
296+
_timeidx = 0
297+
else:
298+
_timeidx = timeidx
299+
348300
if pivot_point is not None:
349301
if pivot_point.lat is not None and pivot_point.lon is not None:
350302
xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx,
@@ -512,30 +464,6 @@ def interpline(field2d, pivot_point=None,
512464
513465
514466
"""
515-
if timeidx is None:
516-
if (latlon is True or is_latlon_pair(start_point) or
517-
is_latlon_pair(pivot_point)):
518-
519-
if wrfin is not None:
520-
# Moving nests aren't supported with ALL_TIMES because the
521-
# domain could move outside of the line, which causes
522-
# crashes or different line lengths.
523-
if is_moving_domain(wrfin):
524-
raise ValueError("Requesting all times with a moving nest "
525-
"is not supported when using a lat/lon "
526-
"line because the domain could "
527-
"move outside of line. "
528-
"You must request each time "
529-
"individually.")
530-
else:
531-
_timeidx = 0
532-
533-
# If using grid coordinates, then don't care about lat/lon
534-
# coordinates. Just use 0.
535-
else:
536-
_timeidx = 0
537-
else:
538-
_timeidx = timeidx
539467

540468
try:
541469
xy = cache["xy"]
@@ -544,6 +472,32 @@ def interpline(field2d, pivot_point=None,
544472
end_point_xy = None
545473
pivot_point_xy = None
546474

475+
if timeidx is None:
476+
if (latlon is True or is_latlon_pair(start_point) or
477+
is_latlon_pair(pivot_point)):
478+
479+
if wrfin is not None:
480+
# Moving nests aren't supported with ALL_TIMES because the
481+
# domain could move outside of the line, which causes
482+
# crashes or different line lengths.
483+
if is_moving_domain(wrfin):
484+
raise ValueError("Requesting all times with a moving nest "
485+
"is not supported when using a lat/lon "
486+
"line because the domain could "
487+
"move outside of line. "
488+
"You must request each time "
489+
"individually.")
490+
else:
491+
# Domain not moving, just use 0
492+
_timeidx = 0
493+
494+
# If using grid coordinates, then don't care about lat/lon
495+
# coordinates. Just use 0.
496+
else:
497+
_timeidx = 0
498+
else:
499+
_timeidx = timeidx
500+
547501
if pivot_point is not None:
548502
if pivot_point.lat is not None and pivot_point.lon is not None:
549503
xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx,
@@ -566,10 +520,10 @@ def interpline(field2d, pivot_point=None,
566520
end_point_xy = (xy_coords.x, xy_coords.y)
567521
else:
568522
end_point_xy = (end_point.x, end_point.y)
569-
523+
570524
xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy,
571525
end_point_xy)
572-
526+
573527
return _interpline(field2d, xy)
574528

575529

‎src/wrf/metadecorators.py‎

Lines changed: 60 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,8 @@
99
from .extension import _interpline
1010
from .util import (extract_vars, either, from_args, arg_location,
1111
is_coordvar, latlon_coordvars, to_np,
12-
from_var, iter_left_indexes, is_mapping)
12+
from_var, iter_left_indexes, is_mapping,
13+
is_moving_domain, is_latlon_pair)
1314
from .coordpair import CoordPair
1415
from .py3compat import viewkeys, viewitems, py3range
1516
from .interputils import get_xy_z_params, get_xy, to_xy_coords
@@ -932,24 +933,49 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
932933
end_point_xy = None
933934
pivot_point_xy = None
934935

936+
if timeidx is None:
937+
if (inc_latlon is True or is_latlon_pair(start_point) or
938+
is_latlon_pair(pivot_point)):
939+
940+
if wrfin is not None:
941+
# Moving nests aren't supported with ALL_TIMES because the
942+
# domain could move outside of the cross section, which causes
943+
# crashes or different line lengths.
944+
if is_moving_domain(wrfin):
945+
raise ValueError("Requesting all times with a moving nest "
946+
"is not supported when using lat/lon "
947+
"cross sections because the domain could "
948+
"move outside of the cross section. "
949+
"You must request each time "
950+
"individually.")
951+
else:
952+
_timeidx = 0
953+
954+
# If using grid coordinates, then don't care about lat/lon
955+
# coordinates. Just use 0.
956+
else:
957+
_timeidx = 0
958+
else:
959+
_timeidx = timeidx
960+
935961
if pivot_point is not None:
936962
if pivot_point.lat is not None and pivot_point.lon is not None:
937-
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
963+
xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx,
938964
stagger, projection, ll_point)
939965
pivot_point_xy = (xy_coords.x, xy_coords.y)
940966
else:
941967
pivot_point_xy = (pivot_point.x, pivot_point.y)
942968

943969
if start_point is not None and end_point is not None:
944970
if start_point.lat is not None and start_point.lon is not None:
945-
xy_coords = to_xy_coords(start_point, wrfin, timeidx,
971+
xy_coords = to_xy_coords(start_point, wrfin, _timeidx,
946972
stagger, projection, ll_point)
947973
start_point_xy = (xy_coords.x, xy_coords.y)
948974
else:
949975
start_point_xy = (start_point.x, start_point.y)
950976

951977
if end_point.lat is not None and end_point.lon is not None:
952-
xy_coords = to_xy_coords(end_point, wrfin, timeidx,
978+
xy_coords = to_xy_coords(end_point, wrfin, _timeidx,
953979
stagger, projection, ll_point)
954980
end_point_xy = (xy_coords.x, xy_coords.y)
955981
else:
@@ -1149,10 +1175,36 @@ def _set_line_meta(wrapped, instance, args, kwargs):
11491175
start_point_xy = None
11501176
end_point_xy = None
11511177
pivot_point_xy = None
1178+
1179+
if timeidx is None:
1180+
if (inc_latlon is True or is_latlon_pair(start_point) or
1181+
is_latlon_pair(pivot_point)):
1182+
1183+
if wrfin is not None:
1184+
# Moving nests aren't supported with ALL_TIMES because the
1185+
# domain could move outside of the line, which causes
1186+
# crashes or different line lengths.
1187+
if is_moving_domain(wrfin):
1188+
raise ValueError("Requesting all times with a moving nest "
1189+
"is not supported when using a lat/lon "
1190+
"line because the domain could "
1191+
"move outside of line. "
1192+
"You must request each time "
1193+
"individually.")
1194+
else:
1195+
# Domain not moving, just use 0
1196+
_timeidx = 0
1197+
1198+
# If using grid coordinates, then don't care about lat/lon
1199+
# coordinates. Just use 0.
1200+
else:
1201+
_timeidx = 0
1202+
else:
1203+
_timeidx = timeidx
11521204

11531205
if pivot_point is not None:
11541206
if pivot_point.lat is not None and pivot_point.lon is not None:
1155-
xy_coords = to_xy_coords(pivot_point, wrfin, timeidx,
1207+
xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx,
11561208
stagger, projection, ll_point)
11571209
pivot_point_xy = (xy_coords.x, xy_coords.y)
11581210
else:
@@ -1161,19 +1213,20 @@ def _set_line_meta(wrapped, instance, args, kwargs):
11611213

11621214
if start_point is not None and end_point is not None:
11631215
if start_point.lat is not None and start_point.lon is not None:
1164-
xy_coords = to_xy_coords(start_point, wrfin, timeidx,
1216+
xy_coords = to_xy_coords(start_point, wrfin, _timeidx,
11651217
stagger, projection, ll_point)
11661218
start_point_xy = (xy_coords.x, xy_coords.y)
11671219
else:
11681220
start_point_xy = (start_point.x, start_point.y)
11691221

11701222
if end_point.lat is not None and end_point.lon is not None:
1171-
xy_coords = to_xy_coords(end_point, wrfin, timeidx,
1223+
xy_coords = to_xy_coords(end_point, wrfin, _timeidx,
11721224
stagger, projection, ll_point)
11731225
end_point_xy = (xy_coords.x, xy_coords.y)
11741226
else:
11751227
end_point_xy = (end_point.x, end_point.y)
11761228

1229+
11771230
xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy, end_point_xy)
11781231

11791232
# Make a copy so we don't modify a user supplied cache

‎test/ncl_get_var.ncl‎

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -236,26 +236,31 @@
236236

237237
t2_line2 = wrf_user_interpline(t2, pivot, xopt)
238238

239+
fout->t2_line2 = t2_line2
240+
239241
lats = wrf_user_getvar(input_file, "lat", 0)
240242
lons = wrf_user_getvar(input_file, "lon", 0)
241243

242244
start_lat = min(lats) + .25d*(max(lats) - min(lats))
243-
end_lat = min(lats) + .75d*(max(lats) - min(lats))
245+
end_lat = min(lats) + .65d*(max(lats) - min(lats))
244246

245247
start_lon = min(lons) + .25d*(max(lons) - min(lons))
246-
end_lon = min(lons) + .75d*(max(lons) - min(lons))
248+
end_lon = min(lons) + .65d*(max(lons) - min(lons))
247249

248250
start_end = (/ start_lon, start_lat, end_lon, end_lat /)
249251

250-
; For the new cross section routine
252+
; For the new line routine
251253
xopt := True
252254
xopt@use_pivot = False
253255
xopt@latlon = True
254256
xopt@file_handle = input_file
255257
xopt@timeidx = 0
256258
xopt@linecoords = True
257-
259+
258260
t2_line3 = wrf_user_interpline(t2, start_end, xopt)
261+
t2_line3!1 = "line_idx_t2_line3"
262+
263+
fout->t2_line3 = t2_line3
259264

260265
times = wrf_user_getvar(input_file, "times", -1)
261266
ntimes = dimsizes(times)
@@ -270,6 +275,22 @@
270275
fout->$name$ = t2_line
271276
end do
272277

278+
; Make sure the 1 time case still works
279+
t2 := wrf_user_getvar(input_file, "T2", 0)
280+
281+
; For the new line routine
282+
xopt := True
283+
xopt@use_pivot = False
284+
xopt@latlon = True
285+
xopt@file_handle = input_file
286+
xopt@timeidx = 0
287+
xopt@linecoords = True
288+
289+
t2_line4 = wrf_user_interpline(t2, start_end, xopt)
290+
t2_line4!0 = "t2_line4_idx"
291+
292+
fout->t2_line4 = t2_line4
293+
273294

274295
;;;;;;;;;;;;;;;;;;;;;;; 3D interpolation to new vertical coordinates
275296
time = -1

0 commit comments

Comments
 (0)