diff --git a/src/wrf/g_latlon.py b/src/wrf/g_latlon.py index 6b1d2da..09e361b 100755 --- a/src/wrf/g_latlon.py +++ b/src/wrf/g_latlon.py @@ -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(): @@ -386,7 +387,128 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0, return _ll_to_xy(latitude, longitude, _wrfin, timeidx, stagger, "cat", 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) @@ -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, *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,11 +603,13 @@ 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, 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, 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, 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) diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 23fb4ca..d875fd1 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -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. diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index e22c1a4..505b204 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -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, diff --git a/test/test_proj_params.py b/test/test_proj_params.py new file mode 100644 index 0000000..77a0ec7 --- /dev/null +++ b/test/test_proj_params.py @@ -0,0 +1,314 @@ +import unittest as ut +import numpy.testing as nt +import numpy as np +import numpy.ma as ma +import os, sys +import subprocess + +from wrf import xy_to_ll_proj, ll_to_xy_proj, to_np + + +class WRFLatLonProjTest(ut.TestCase): + longMessage = True + +def make_test(xy_or_ll_out): + def test(self): + + # Python 3 + if sys.version_info > (3,): + assert_raises_regex = self.assertRaisesRegex + xrange = range + else: + assert_raises_regex = self.assertRaisesRegexp + + if xy_or_ll_out == "xy": + # Test the required failures + with assert_raises_regex(ValueError, ".*map_proj.*" ): + ll_to_xy_proj(30,-110) + + with assert_raises_regex(ValueError, ".*ref_lat.*" ): + ll_to_xy_proj(30,-110, map_proj=1) + + with assert_raises_regex(ValueError, ".*ref_lon.*" ): + ll_to_xy_proj(30,-110, map_proj=1, ref_lat=45) + + with assert_raises_regex(ValueError, ".*known_x.*" ): + ll_to_xy_proj(30,-110, map_proj=1, ref_lat=45.0, + ref_lon=-120.) + + with assert_raises_regex(ValueError, ".*known_y.*" ): + ll_to_xy_proj(30,-110, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=1) + + with assert_raises_regex(ValueError, ".*dx.*" ): + ll_to_xy_proj(30,-110, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0) + + ####### Now test the projections + + # Lambert Conformal - truelat1, stand_lon required + with assert_raises_regex(ValueError, ".*truelat1.*" ): + ll_to_xy_proj(30,-110, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000.) + + with assert_raises_regex(ValueError, ".*stand_lon.*" ): + ll_to_xy_proj(30,-110, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000., + truelat1=60.) + + # Make sure it runs with all params set vs not + p_all = ll_to_xy_proj(28., -89., map_proj=1, ref_lat=17.803, + ref_lon=-100.7747, known_x=0, known_y=0, + dx=30000., truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=90., pole_lon=0.) + + + p_def = ll_to_xy_proj(28., -89., map_proj=1, ref_lat=17.803, + ref_lon=-100.7747, known_x=0, known_y=0, + dx=30000., truelat1=30., + stand_lon=-89.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + # Polar Stereographic - truelat1, stand_lon + + with assert_raises_regex(ValueError, ".*truelat1.*" ): + ll_to_xy_proj(30,-110, map_proj=2, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000.) + + with assert_raises_regex(ValueError, ".*stand_lon.*" ): + ll_to_xy_proj(30,-110, map_proj=2, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000., + truelat1=60.) + + p_all = ll_to_xy_proj(28., -89., map_proj=2, ref_lat=17.933, + ref_lon=-100.0735, known_x=0, known_y=0, + dx=30000., truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=90., pole_lon=0.) + + + p_def = ll_to_xy_proj(28., -89., map_proj=2, ref_lat=17.933, + ref_lon=-100.0735, known_x=0, known_y=0, + dx=30000., truelat1=30., + stand_lon=-89.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + + # Mercator - truelat1 + + with assert_raises_regex(ValueError, ".*truelat1.*" ): + ll_to_xy_proj(30,-110, map_proj=2, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000.) + + p_all = ll_to_xy_proj(28., -89., map_proj=3, ref_lat=19.1075, + ref_lon=-101.008, known_x=0, known_y=0, + dx=30000., truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=90., pole_lon=0.) + + + p_def = ll_to_xy_proj(28., -89., map_proj=3, ref_lat=19.1075, + ref_lon=-101.008, known_x=0, known_y=0, + dx=30000., truelat1=30.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + # Lat/lon - stand_lon, dy, pole_lat, pole_lon + + with assert_raises_regex(ValueError, ".*stand_lon.*" ): + ll_to_xy_proj(30,-110, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=.2698388) + + with assert_raises_regex(ValueError, ".*dy.*" ): + ll_to_xy_proj(30,-110, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + stand_lon=89.0, + dx=.2698388) + + with assert_raises_regex(ValueError, ".*pole_lat.*" ): + ll_to_xy_proj(30,-110, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + stand_lon=89.0, + dx=.2698388, dy=.2698388) + + with assert_raises_regex(ValueError, ".*pole_lon.*" ): + ll_to_xy_proj(30,-110, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + stand_lon=89.0, + dx=.2698388, dy=.2698388, + pole_lat=62.0) + + p_all = ll_to_xy_proj(28.,-89., map_proj=6, ref_lat=17.6759, + ref_lon=-101.4286, known_x=0, known_y=0, + dx=30000, dy=30000, + truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=62.0, + pole_lon=180.) + + + p_def = ll_to_xy_proj(28., -89., map_proj=6, ref_lat=17.6759, + ref_lon=-101.4286, known_x=0, known_y=0, + stand_lon=-89., + dx=30000, dy=30000, pole_lat=62.0, + pole_lon=180.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + + if xy_or_ll_out == "ll": + + # Test the required failures + with assert_raises_regex(ValueError, ".*map_proj.*" ): + xy_to_ll_proj(45, 50) + + with assert_raises_regex(ValueError, ".*ref_lat.*" ): + xy_to_ll_proj(45, 50, map_proj=1) + + with assert_raises_regex(ValueError, ".*ref_lon.*" ): + xy_to_ll_proj(45, 50, map_proj=1, ref_lat=45) + + with assert_raises_regex(ValueError, ".*known_x.*" ): + xy_to_ll_proj(45, 50, map_proj=1, ref_lat=45.0, + ref_lon=-120.) + + with assert_raises_regex(ValueError, ".*known_y.*" ): + xy_to_ll_proj(45, 50, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=1) + + with assert_raises_regex(ValueError, ".*dx.*" ): + xy_to_ll_proj(45, 50, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0) + + ####### Now test the projections + + # Lambert Conformal - truelat1, stand_lon required + with assert_raises_regex(ValueError, ".*truelat1.*" ): + xy_to_ll_proj(45, 50, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000.) + + with assert_raises_regex(ValueError, ".*stand_lon.*" ): + xy_to_ll_proj(45, 50, map_proj=1, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000., + truelat1=60.) + + # Make sure it runs with all params set vs not + p_all = xy_to_ll_proj(45, 50, map_proj=1, ref_lat=17.803, + ref_lon=-100.7747, known_x=0, known_y=0, + dx=30000., truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=90., pole_lon=0.) + + + p_def = xy_to_ll_proj(45, 50, map_proj=1, ref_lat=17.803, + ref_lon=-100.7747, known_x=0, known_y=0, + dx=30000., truelat1=30., + stand_lon=-89.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + # Polar Stereographic - truelat1, stand_lon + + with assert_raises_regex(ValueError, ".*truelat1.*" ): + xy_to_ll_proj(45, 50, map_proj=2, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000.) + + with assert_raises_regex(ValueError, ".*stand_lon.*" ): + xy_to_ll_proj(45, 50, map_proj=2, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000., + truelat1=60.) + + p_all = xy_to_ll_proj(45, 50, map_proj=2, ref_lat=17.933, + ref_lon=-100.0735, known_x=0, known_y=0, + dx=30000., truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=90., pole_lon=0.) + + + p_def = xy_to_ll_proj(45, 50, map_proj=2, ref_lat=17.933, + ref_lon=-100.0735, known_x=0, known_y=0, + dx=30000., truelat1=30., + stand_lon=-89.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + + # Mercator - truelat1 + + with assert_raises_regex(ValueError, ".*truelat1.*" ): + xy_to_ll_proj(45, 50, map_proj=2, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=3000.) + + p_all = xy_to_ll_proj(45, 50, map_proj=3, ref_lat=19.1075, + ref_lon=-101.008, known_x=0, known_y=0, + dx=30000., truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=90., pole_lon=0.) + + + p_def = xy_to_ll_proj(45, 50, map_proj=3, ref_lat=19.1075, + ref_lon=-101.008, known_x=0, known_y=0, + dx=30000., truelat1=30.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + # Lat/lon - stand_lon, dy, pole_lat, pole_lon + + with assert_raises_regex(ValueError, ".*stand_lon.*" ): + xy_to_ll_proj(45, 50, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + dx=.2698388) + + with assert_raises_regex(ValueError, ".*dy.*" ): + xy_to_ll_proj(45, 50, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + stand_lon=89.0, + dx=.2698388) + + with assert_raises_regex(ValueError, ".*pole_lat.*" ): + xy_to_ll_proj(45, 50, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + stand_lon=89.0, + dx=.2698388, dy=.2698388) + + with assert_raises_regex(ValueError, ".*pole_lon.*" ): + xy_to_ll_proj(45, 50, map_proj=6, ref_lat=45.0, + ref_lon=-120., known_x=0, known_y=0, + stand_lon=89.0, + dx=.2698388, dy=.2698388, + pole_lat=62.0) + + p_all = xy_to_ll_proj(64, 40, map_proj=6, ref_lat=17.6759, + ref_lon=-101.4286, known_x=0, known_y=0, + dx=30000, dy=30000, + truelat1=30., truelat2=30., + stand_lon=-89., pole_lat=62.0, + pole_lon=180.) + + + p_def = xy_to_ll_proj(64, 40, map_proj=6, ref_lat=17.6759, + ref_lon=-101.4286, known_x=0, known_y=0, + stand_lon=-89., + dx=30000, dy=30000, pole_lat=62.0, + pole_lon=180.) + + nt.assert_allclose(to_np(p_all), to_np(p_def)) + + return test + + +if __name__ == "__main__": + + for v in ("xy", "ll"): + test_func = make_test(v) + setattr(WRFLatLonProjTest, 'test_{0}'.format(v), test_func) + + ut.main() + \ No newline at end of file diff --git a/test/utests.py b/test/utests.py index 613b13e..114d1a4 100644 --- a/test/utests.py +++ b/test/utests.py @@ -47,7 +47,7 @@ def setUpModule(): PATTERN, outfile) - print cmd + print(cmd) if not os.path.exists(outfile): status = subprocess.call(cmd, shell=True) @@ -179,8 +179,8 @@ def make_test(varname, dir, pattern, referent, multi=False, pynio=False): except: absdiff = np.abs(to_np(my_vals) - ref_vals) maxdiff = np.amax(absdiff) - print (maxdiff) - print np.argwhere(absdiff == maxdiff) + print(maxdiff) + print(np.argwhere(absdiff == maxdiff)) raise @@ -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.) @@ -890,7 +892,7 @@ def make_latlon_test(testid, dir, pattern, referent, single, y_s = np.asarray([10, 50, 90], int) if single: - timeidx=8 + timeidx = 8 ref_vals = refnc.variables["ll2"][:] ll = xy_to_ll(wrfin, x_s[0], y_s[0], timeidx=timeidx) ref = ref_vals[::-1,0] @@ -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)