Browse Source

ll_to_xy_proj and xy_to_ll_proj now check for required args and set defaults.

Fixes #61.
lon0
Bill Ladwig 7 years ago
parent
commit
d15d4a3e2d
  1. 142
      src/wrf/g_latlon.py
  2. 5
      src/wrf/interp.py
  3. 1
      src/wrf/latlonutils.py
  4. 10
      test/utests.py

142
src/wrf/g_latlon.py

@ -8,6 +8,7 @@ from .util import extract_vars, get_id, get_iterable, is_mapping, to_np @@ -8,6 +8,7 @@ from .util import extract_vars, get_id, get_iterable, is_mapping, to_np
from .py3compat import viewkeys
from .latlonutils import _lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll
from .metadecorators import set_latlon_metadata
from .constants import Constants, ProjectionTypes
from .config import xarray_enabled
if xarray_enabled():
@ -388,6 +389,127 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0, @@ -388,6 +389,127 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0,
squeeze, None, _key, as_int, **{})
def _set_defaults(projparams):
"""Check projection parameters and set defaults.
Throws an exception if projection parameters required by WPS are not
provided, along with any other parameters required by the map projection
transformation routines.
For parameters not used by the projection type, defaults are used so that
the None values don't pass through to fortran downstream.
Args:
projparams (:obj:`dict`): Projection parameters dictionary.
Returns:
:obj:`dict`: The projection parameters with default values inserted
where applicable.
"""
_params = dict(projparams)
map_proj = _params.get("map_proj")
# All projections require these arguments
if map_proj is None:
raise ValueError("'map_proj' cannot be None")
if _params.get("ref_lat") is None:
raise ValueError("'ref_lat' cannot be None")
if _params.get("ref_lon") is None:
raise ValueError("'ref_lon' cannot be None")
if _params.get("known_x") is None:
raise ValueError("'known_x' cannot be None")
if _params.get("known_y") is None:
raise ValueError("'known_y' cannot be None")
if _params.get("dx") is None:
raise ValueError("'dx' cannot be None")
# Requires truelat1,stand_lon, truelat2, dx, dy
if map_proj == ProjectionTypes.LAMBERT_CONFORMAL:
if _params.get("truelat1") is None:
raise ValueError("'truelat1' cannot be None")
if _params.get("stand_lon") is None:
raise ValueError("'stand_lon' cannot be None")
if _params.get("truelat2") is None:
_params["truelat2"] = _params["truelat1"]
# Requires truelat1, stand_lon
elif map_proj == ProjectionTypes.POLAR_STEREOGRAPHIC:
if _params.get("truelat1") is None:
raise ValueError("'truelat1' cannot be None")
if _params.get("stand_lon") is None:
raise ValueError("'stand_lon' cannot be None")
# Requires truelat1
elif map_proj == ProjectionTypes.MERCATOR:
if _params.get("truelat1") is None:
raise ValueError("'truelat1' cannot be None")
if _params.get("stand_lon") is None:
_params["stand_lon"] = 0.0
# Requires pole_lat, pole_lon, stand_lon
elif map_proj == ProjectionTypes.LAT_LON:
if _params.get("stand_lon") is None:
raise ValueError("'stand_lon' cannot be None")
if _params.get("dy") is None:
raise ValueError("'dy' cannot be None")
if _params.get("pole_lat") is None:
raise ValueError("'pole_lat' cannot be None")
if _params.get("pole_lon") is None:
raise ValueError("'pole_lon' cannot be None")
if _params.get("latinc") is None:
_params["latinc"] = ((_params["dy"]*360.0)/2.0 /
Constants.PI/Constants.WRF_EARTH_RADIUS)
if _params.get("loninc") is None:
_params["loninc"] = ((_params["dx"]*360.0)/2.0 /
Constants.PI/Constants.WRF_EARTH_RADIUS)
else:
raise ValueError("invalid 'map_proj' value of {}".format(map_proj))
# Set these to defaults if not used so that the Fortran routines
# don't crash
if _params.get("truelat1") is None:
_params["truelat1"] = 0.
if _params.get("truelat2") is None:
_params["truelat2"] = 0.
if _params.get("pole_lat") is None:
_params["pole_lat"] = 90.0
if _params.get("pole_lon") is None:
_params["pole_lon"] = 0.0
if _params.get("dx") is None:
_params["dx"] = 0.0
if _params.get("dy") is None:
_params["dy"] = 0.0
if _params.get("latinc") is None:
_params["latinc"] = 0.
if _params.get("loninc") is None:
_params["loninc"] = 0.
return _params
@set_latlon_metadata(xy=True)
def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
@ -435,7 +557,8 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True, @@ -435,7 +557,8 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
truelat2 (:obj:`float`): True latitude 2. Optional for
map_proj = 1 (defaults to 0 otherwise).
stand_lon (:obj:`float`): Standard longitude. Required.
stand_lon (:obj:`float`): Standard longitude. Required for map_proj =
1, 2, 6 (defaults to 0 otherwise).
ref_lat (:obj:`float`): A reference latitude. Required.
@ -454,10 +577,10 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True, @@ -454,10 +577,10 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
*map_proj* = 6 (defaults to 0 otherwise).
dx (:obj:`float`): The x spacing in meters at the true latitude.
Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise).
Required for all map projections.
dy (:obj:`float`) - The y spacing in meters at the true latitude.
Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise).
Required for *map_proj* = 6 (defaults to 0 otherwise).
latinc (:obj:`float`): Required for *map_proj* = 6. Defined as:
@ -480,12 +603,14 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True, @@ -480,12 +603,14 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
"""
loc = locals()
projparams = {name : loc[name] for name in ("map_proj", "truelat1",
_projparams = {name : loc[name] for name in ("map_proj", "truelat1",
"truelat2", "stand_lon", "ref_lat",
"ref_lon", "pole_lat", "pole_lon",
"known_x", "known_y", "dx", "dy",
"latinc", "loninc")}
projparams = _set_defaults(_projparams)
return _ll_to_xy(latitude, longitude, None, 0, True, "cat", squeeze, None,
None, as_int, **projparams)
@ -600,7 +725,8 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None, @@ -600,7 +725,8 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None,
truelat2 (:obj:`float`): True latitude 2. Optional for
map_proj = 1 (defaults to 0 otherwise).
stand_lon (:obj:`float`): Standard longitude. Required.
stand_lon (:obj:`float`): Standard longitude. Required for map_proj =
1, 2, 6 (defaults to 0 otherwise).
ref_lat (:obj:`float`): A reference latitude. Required.
@ -622,7 +748,7 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None, @@ -622,7 +748,7 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None,
Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise).
dy (:obj:`float`) - The y spacing in meters at the true latitude.
Required for *map_proj* = 1, 2, 3 (defaults to 0 otherwise).
Required for *map_proj* = 6 (defaults to 0 otherwise).
latinc (:obj:`float`): Required for *map_proj* = 6. Defined as:
@ -645,11 +771,13 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None, @@ -645,11 +771,13 @@ def xy_to_ll_proj(x, y, meta=True, squeeze=True, map_proj=None, truelat1=None,
be a :class:`numpy.ndarray` object with no metadata.
"""
loc = locals()
projparams = {name : loc[name] for name in ("map_proj", "truelat1",
_projparams = {name : loc[name] for name in ("map_proj", "truelat1",
"truelat2", "stand_lon", "ref_lat",
"ref_lon", "pole_lat", "pole_lon",
"known_x", "known_y", "dx", "dy",
"latinc", "loninc")}
projparams = _set_defaults(_projparams)
return _xy_to_ll(x, y, None, 0, None, "cat", squeeze, None, None,
**projparams)

5
src/wrf/interp.py

@ -350,10 +350,9 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), @@ -350,10 +350,9 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
@set_interp_metadata("line")
def interpline(field2d, pivot_point=None,
wrfin=None, timeidx=0, stagger=None, projection=None,
def interpline(field2d, wrfin=None, timeidx=0, stagger=None, projection=None,
ll_point=None,
angle=None, start_point=None,
pivot_point=None, angle=None, start_point=None,
end_point=None, latlon=False,
cache=None, meta=True):
"""Return the two-dimensional field interpolated along a line.

1
src/wrf/latlonutils.py

@ -233,7 +233,6 @@ def _kwarg_proj_params(**projparams): @@ -233,7 +233,6 @@ def _kwarg_proj_params(**projparams):
# Sanity checks
# Required args for all projections
for name, var in viewitems({"MAP_PROJ" : map_proj,
"STAND_LON" : stdlon,
"REF_LAT" : ref_lat,
"REF_LON" : ref_lon,
"KNOWN_X" : known_x,

10
test/utests.py

@ -47,7 +47,7 @@ def setUpModule(): @@ -47,7 +47,7 @@ def setUpModule():
PATTERN,
outfile)
print cmd
print(cmd)
if not os.path.exists(outfile):
status = subprocess.call(cmd, shell=True)
@ -180,7 +180,7 @@ def make_test(varname, dir, pattern, referent, multi=False, pynio=False): @@ -180,7 +180,7 @@ def make_test(varname, dir, pattern, referent, multi=False, pynio=False):
absdiff = np.abs(to_np(my_vals) - ref_vals)
maxdiff = np.amax(absdiff)
print(maxdiff)
print np.argwhere(absdiff == maxdiff)
print(np.argwhere(absdiff == maxdiff))
raise
@ -519,8 +519,10 @@ def make_interp_test(varname, dir, pattern, referent, multi=False, @@ -519,8 +519,10 @@ def make_interp_test(varname, dir, pattern, referent, multi=False,
int(lats.shape[-1]/2)],
lon=lons[int(lons.shape[-2]/2),
int(lons.shape[-1]/2)])
l1 = interpline(t2,wrfin,pivot_point=pivot_point,
l1 = interpline(t2,wrfin=wrfin,pivot_point=pivot_point,
angle=90.0)
l2 = interpline(t2,projection=t2.attrs["projection"],
ll_point=ll_point,
pivot_point=pivot_point, angle=90.)
@ -902,6 +904,8 @@ def make_latlon_test(testid, dir, pattern, referent, single, @@ -902,6 +904,8 @@ def make_latlon_test(testid, dir, pattern, referent, single,
ll_proj = xy_to_ll_proj(x_s[0], y_s[0], **projparams)
nt.assert_allclose(to_np(ll_proj), to_np(ll))
else:
ref_vals = refnc.variables["ll1"][:]
ll = xy_to_ll(wrfin, x_s, y_s, timeidx=None)

Loading…
Cancel
Save