Browse Source

Merge branch 'bugfix/issue_61' into develop

lon0
Bill Ladwig 6 years ago
parent
commit
243158ec81
  1. 142
      src/wrf/g_latlon.py
  2. 5
      src/wrf/interp.py
  3. 1
      src/wrf/latlonutils.py
  4. 314
      test/test_proj_params.py
  5. 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,

314
test/test_proj_params.py

@ -0,0 +1,314 @@ @@ -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()

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