From b8bc1fbb384c65f5b0a66da223d517b29d1d3904 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Tue, 23 Jan 2018 16:29:27 -0700 Subject: [PATCH 1/2] Fixes vertcross and interpline when using lat/lon coordinates with a projection object. The fix requires the user to supply the lower left lat/lon corner points for the domain. Added unit tests. Updated docs. Fixes #44. --- doc/source/basic_usage.rst | 12 ++++++++---- src/wrf/interp.py | 30 ++++++++++++++++++++++++++---- src/wrf/interputils.py | 24 +++++++++++++++++++----- src/wrf/metadecorators.py | 18 ++++++++++++------ test/utests.py | 32 +++++++++++++++++++++++++++++++- 5 files changed, 96 insertions(+), 20 deletions(-) diff --git a/doc/source/basic_usage.rst b/doc/source/basic_usage.rst index 145eed4..c29f03f 100644 --- a/doc/source/basic_usage.rst +++ b/doc/source/basic_usage.rst @@ -736,8 +736,9 @@ Example Using Lat/Lon Coordinates start_point = CoordPair(lat=start_lat, lon=start_lon) end_point = CoordPair(lat=end_lat, lon=end_lon) - # When using lat/lon coordinates, you must supply a netcdf file object, or a - # projection object. + # When using lat/lon coordinates, you must supply a WRF netcdf file object, + # or a projection object with the lower left latitude and lower left + # longitude points. p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True) print(p_vert) @@ -983,8 +984,11 @@ Example Using Lat/Lon Coordinates start_point = CoordPair(lat=start_lat, lon=start_lon) end_point = CoordPair(lat=end_lat, lon=end_lon) - # Calculate the vertical cross section. By setting latlon to True, this - # also calculates the latitude and longitude coordinates along the line + # Calculate the interpolated line. To use latitude and longitude points, + # you must supply a WRF NetCDF file object, or a projection object along + # with the lower left latitude and lower left longitude points. + # Also, by setting latlon to True, this routine will + # also calculate the latitude and longitude coordinates along the line # and adds them to the metadata to help with plotting labels. t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True) diff --git a/src/wrf/interp.py b/src/wrf/interp.py index fcd3a8e..9e5ead6 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -91,7 +91,8 @@ def interplevel(field3d, vert, desiredlev, missing=default_fill(np.float64), @set_interp_metadata("cross") def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), - wrfin=None, timeidx=0, stagger=None, projection=None, + wrfin=None, timeidx=0, stagger=None, projection=None, + ll_lat=None, ll_lon=None, pivot_point=None, angle=None, start_point=None, end_point=None, latlon=False, cache=None, meta=True): @@ -165,6 +166,16 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64), coordinates, and must be specified if *wrfin* is None. Default is None. + ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower + left latitude(s) for your domain, and must be specified if + *wrfin* is None. If the domain is a moving nest, this should be a + sequence of lower left latitudes. Default is None. + + ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower + left longitude(s) for your domain, and must be specified if + *wrfin* is None. If the domain is a moving nest, this should be + a sequence of lower left longitudes. Default is None. + pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for the pivot point, which indicates the location through which the plane will pass. Must also specify *angle*. The coordinate @@ -291,6 +302,7 @@ 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, + ll_lat=None, ll_lon=None, angle=None, start_point=None, end_point=None, latlon=False, cache=None, meta=True): @@ -340,6 +352,16 @@ def interpline(field2d, pivot_point=None, not be used when working with x,y coordinates. Default is None. + ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower + left latitude(s) for your domain, and must be specified if + *wrfin* is None. If the domain is a moving nest, this should be a + sequence of lower left latitudes. Default is None. + + ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower + left longitude(s) for your domain, and must be specified if + *wrfin* is None. If the domain is a moving nest, this should be + a sequence of lower left longitudes. Default is None. + pivot_point (:class:`wrf.CoordPair`, optional): A coordinate pair for the pivot point, which indicates the location through which the plane will pass. Must also specify *angle*. The coordinate @@ -420,7 +442,7 @@ def interpline(field2d, pivot_point=None, if pivot_point is not None: if pivot_point.lat is not None and pivot_point.lon is not None: xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -428,14 +450,14 @@ def interpline(field2d, pivot_point=None, if start_point is not None and end_point is not None: if start_point.lat is not None and start_point.lon is not None: xy_coords = to_xy_coords(start_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) start_point_xy = (xy_coords.x, xy_coords.y) else: start_point_xy = (start_point.x, start_point.y) if end_point.lat is not None and end_point.lon is not None: xy_coords = to_xy_coords(end_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index fb2768f..5289f08 100644 --- a/src/wrf/interputils.py +++ b/src/wrf/interputils.py @@ -329,7 +329,8 @@ def get_xy(var, pivot_point=None, angle=None, return xy -def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None): +def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None, + ll_lat=None, ll_lon=None): """Return the coordinate pairs in grid space. This function converts latitude,longitude coordinate pairs to @@ -371,6 +372,16 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None): projection object to use when working with latitude, longitude coordinates, and must be specified if *wrfin* is None. Default is None. + + ll_lat (:obj:`float`, sequence of :obj:`float`, optional): The lower + left latitude(s) for your domain, and must be specified if + *wrfin* is None. If the domain is a moving nest, this should be a + sequence of lower left latitudes. Default is None. + + ll_lon (:obj:`float`, sequence of :obj:`float`, optional): The lower + left longitude(s) for your domain, and must be specified if + *wrfin* is None. If the domain is a moving nest, this should be + a sequence of lower left longitudes. Default is None. Returns: @@ -379,8 +390,11 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None): """ - if wrfin is None and projection is None: - raise ValueError ("'wrfin' or 'projection' parameter is required") + if (wrfin is None and + (projection is None or ll_lat is None or ll_lon is None)): + raise ValueError ("'wrfin' parameter or " + "'projection', 'll_lat', and 'll_lon' parameters " + "are required") if isinstance(pairs, Iterable): lat = [pair.lat for pair in pairs] @@ -415,8 +429,8 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None): truelat1=projection.truelat1, truelat2=projection.truelat2, stand_lon=projection.stand_lon, - ref_lat=projection.ll_lat, - ref_lon=projection.ll_lon, + ref_lat=ll_lat, + ref_lon=ll_lon, pole_lat=pole_lat, pole_lon=pole_lon, known_x=0, diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 228c857..2abcb63 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -887,6 +887,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): argvars = from_args(wrapped, ("field3d", "vert", "levels", "latlon", "missing", "wrfin", "timeidx", "stagger", "projection", + "ll_lat", "ll_lon", "pivot_point", "angle", "start_point", "end_point", "cache"), @@ -901,6 +902,8 @@ def _set_cross_meta(wrapped, instance, args, kwargs): timeidx = argvars["timeidx"] stagger = argvars["stagger"] projection = argvars["projection"] + ll_lat = argvars["ll_lat"] + ll_lon = argvars["ll_lon"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] start_point = argvars["start_point"] @@ -914,7 +917,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): if pivot_point is not None: if pivot_point.lat is not None and pivot_point.lon is not None: xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -922,14 +925,14 @@ def _set_cross_meta(wrapped, instance, args, kwargs): if start_point is not None and end_point is not None: if start_point.lat is not None and start_point.lon is not None: xy_coords = to_xy_coords(start_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) start_point_xy = (xy_coords.x, xy_coords.y) else: start_point_xy = (start_point.x, start_point.y) if end_point.lat is not None and end_point.lon is not None: xy_coords = to_xy_coords(end_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) @@ -1104,6 +1107,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): """ argvars = from_args(wrapped, ("field2d", "wrfin", "timeidx", "stagger", "projection", + "ll_lat", "ll_lon", "pivot_point", "angle", "start_point", "end_point", "latlon", "cache"), @@ -1114,6 +1118,8 @@ def _set_line_meta(wrapped, instance, args, kwargs): timeidx = argvars["timeidx"] stagger = argvars["stagger"] projection = argvars["projection"] + ll_lat = argvars["ll_lat"] + ll_lon = argvars["ll_lon"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] start_point = argvars["start_point"] @@ -1131,7 +1137,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): if pivot_point is not None: if pivot_point.lat is not None and pivot_point.lon is not None: xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) pivot_point_xy = (xy_coords.x, xy_coords.y) else: pivot_point_xy = (pivot_point.x, pivot_point.y) @@ -1140,14 +1146,14 @@ def _set_line_meta(wrapped, instance, args, kwargs): if start_point is not None and end_point is not None: if start_point.lat is not None and start_point.lon is not None: xy_coords = to_xy_coords(start_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) start_point_xy = (xy_coords.x, xy_coords.y) else: start_point_xy = (start_point.x, start_point.y) if end_point.lat is not None and end_point.lon is not None: xy_coords = to_xy_coords(end_point, wrfin, timeidx, - stagger, projection) + stagger, projection, ll_lat, ll_lon) end_point_xy = (xy_coords.x, xy_coords.y) else: end_point_xy = (end_point.x, end_point.y) diff --git a/test/utests.py b/test/utests.py index 75cf254..d00fef8 100644 --- a/test/utests.py +++ b/test/utests.py @@ -293,7 +293,23 @@ def make_interp_test(varname, wrf_in, referent, multi=False, # Note: Until the bug is fixed in NCL, the wrf-python cross # sections will contain one extra point - nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross, rtol=.01) + nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross, + rtol=.01) + + # Test the manual projection method with lat/lon + lats = hts.coords["XLAT"] + lons = hts.coords["XLONG"] + ll_lat = lats[0,0] + ll_lon = lons[0,0] + pivot = CoordPair(lat=lats[lats.shape[-2]/2, lats.shape[-1]/2], + lon=lons[lons.shape[-2]/2, lons.shape[-1]/2]) + v1 = vertcross(hts,p,wrfin=in_wrfnc,pivot_point=pivot_point, + angle=90.0) + v2 = vertcross(hts,p,projection=hts.attrs["projection"], + ll_lat=ll_lat, ll_lon=ll_lon, + pivot_point=pivot_point, angle=90.) + + nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01) # Test opposite p_cross1 = vertcross(p,hts,pivot_point=pivot_point, angle=90.0) @@ -329,6 +345,20 @@ def make_interp_test(varname, wrf_in, referent, multi=False, # Note: After NCL is fixed, remove the slice. nt.assert_allclose(to_np(t2_line1)[...,0:-1], ref_t2_line) + # Test the manual projection method with lat/lon + lats = t2.coords["XLAT"] + lons = t2.coords["XLONG"] + ll_lat = lats[0,0] + ll_lon = lons[0,0] + pivot = CoordPair(lat=lats[lats.shape[-2]/2, lats.shape[-1]/2], + lon=lons[lons.shape[-2]/2, lons.shape[-1]/2]) + l1 = interpline(t2,wrfin=in_wrfnc,pivot_point=pivot_point, + angle=90.0) + l2 = interpline(t2,projection=t2.attrs["projection"], + ll_lat=ll_lat, ll_lon=ll_lon, + pivot_point=pivot_point, angle=90.) + nt.assert_allclose(to_np(l1), to_np(l2), rtol=.01) + # Test point to point start_point = CoordPair(0, t2.shape[-2]/2) end_point = CoordPair(-1, t2.shape[-2]/2) From 37052d323590170e621eb28d8bdeb01cd027a7a4 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Tue, 23 Jan 2018 16:37:44 -0700 Subject: [PATCH 2/2] Updated documentation --- doc/source/new.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/doc/source/new.rst b/doc/source/new.rst index e2ee510..e2ff9cc 100644 --- a/doc/source/new.rst +++ b/doc/source/new.rst @@ -44,6 +44,12 @@ v1.1.0 vertically staggered grid. - A DOI is now available for wrf-python. Please cite wrf-python if you are using it for your research. (See :ref:`citation`) +- Fixed issue with vertcross and interpline not working correctly when a + projection object is used. Users will now have to supply the lower left + latitude and longitude corner point. +- Beginning with numpy 1.14, wrf-python can be built using the MSVC + compiler with gfortran. WRF-Python can now be built for Python 3.5+ on + services like AppVeyor. v1.0.5