From d1a4651c2346728102ff7788fbbb9ab71ebea80d Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 1 Feb 2019 12:18:45 -0700 Subject: [PATCH] PEP 8. Also reorganized the test directory. Various snippets that may or may not be of use anymore added the the misc directory. Misc. NCL scripts used during development added to ncl directory. --- doc/source/_templates/product_table.txt | 10 + doc/source/internal_api/index.rst | 3 + src/wrf/g_geoht.py | 68 ++++ src/wrf/routines.py | 5 +- test/ctt_test.py | 17 +- test/generator_test.py | 10 +- test/misc/extract_one_time.py | 40 +++ test/misc/loop_and_fill_meta.py | 70 ++++ test/misc/mocktest.py | 43 ++- test/misc/projtest.py | 188 +++++----- test/misc/quiver_test.py | 58 +++ test/misc/reduce_file.py | 14 +- test/misc/snippet.py | 39 +- test/{ => misc}/varcache.py | 34 +- test/{ => misc}/viewtest.py | 19 +- test/misc/wps.py | 230 ++++++++++++ test/{ => ncl}/listBug.ncl | 0 test/{ => ncl}/ncl_get_var.ncl | 0 test/ncl/ncl_vertcross.ncl | 92 +++++ test/ncl/refl10_cross.ncl | 81 +++++ test/ncl/rotated_test.ncl | 26 ++ test/ncl/test_this.ncl | 21 ++ test/{ => ncl}/test_vinterp.ncl | 0 test/ncl/wrf_user_vertcross.ncl | 416 ++++++++++++++++++++++ test/test_filevars.py | 80 +++-- test/test_inputs.py | 103 +++--- test/test_multi_cache.py | 42 ++- test/test_omp.py | 91 +++-- test/test_proj_params.py | 454 ++++++++++++------------ test/test_units.py | 33 +- test/utests.py | 2 +- 31 files changed, 1709 insertions(+), 580 deletions(-) create mode 100644 test/misc/extract_one_time.py create mode 100644 test/misc/loop_and_fill_meta.py create mode 100644 test/misc/quiver_test.py rename test/{ => misc}/varcache.py (68%) rename test/{ => misc}/viewtest.py (51%) create mode 100644 test/misc/wps.py rename test/{ => ncl}/listBug.ncl (100%) rename test/{ => ncl}/ncl_get_var.ncl (100%) create mode 100644 test/ncl/ncl_vertcross.ncl create mode 100644 test/ncl/refl10_cross.ncl create mode 100644 test/ncl/rotated_test.ncl create mode 100644 test/ncl/test_this.ncl rename test/{ => ncl}/test_vinterp.ncl (100%) create mode 100644 test/ncl/wrf_user_vertcross.ncl diff --git a/doc/source/_templates/product_table.txt b/doc/source/_templates/product_table.txt index e3fe65b..2ecc5f1 100644 --- a/doc/source/_templates/product_table.txt +++ b/doc/source/_templates/product_table.txt @@ -245,6 +245,16 @@ | | | | | | | | mi | | +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ +| height_agl | Model Height for Mass Grid (AGL) | m | **units** (str) : Set to desired units. Default is *'m'*. | +| | | | | +| | | km | | +| | | | | +| | | dm | | +| | | | | +| | | ft | | +| | | | | +| | | mi | | ++--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ | zstag | Model Height for Vertically Staggered Grid | m | **msl** (boolean): Set to False to return AGL values. True is for MSL. Default is *True*. | | | | | | | | | km | **units** (str) : Set to desired units. Default is *'m'*. | diff --git a/doc/source/internal_api/index.rst b/doc/source/internal_api/index.rst index cfd3c7c..089bbb0 100644 --- a/doc/source/internal_api/index.rst +++ b/doc/source/internal_api/index.rst @@ -23,6 +23,9 @@ The routines below are called internally by :meth:`wrf.getvar`. wrf.g_dewpoint.get_dp_2m wrf.g_geoht.get_geopt wrf.g_geoht.get_height + wrf.g_geoht.get_height_agl + wrf.g_geoht.get_stag_geopt + wrf.g_geoht.get_stag_height wrf.g_helicity.get_srh wrf.g_helicity.get_uh wrf.g_omega.get_omega diff --git a/src/wrf/g_geoht.py b/src/wrf/g_geoht.py index 7e78e5a..eb0b841 100755 --- a/src/wrf/g_geoht.py +++ b/src/wrf/g_geoht.py @@ -391,3 +391,71 @@ def get_stag_height(wrfin, timeidx=0, method="cat", squeeze=True, return _get_geoht(wrfin, timeidx, method, squeeze, cache, meta, _key, True, msl, stag=True) + + +@set_height_metadata(geopt=False, stag=False) +@convert_units("height", "m") +def get_height_agl(wrfin, timeidx=0, method="cat", squeeze=True, + cache=None, meta=True, _key=None, units="m"): + """Return the geopotential height (AGL). + + The geopotential height is returned as Above Ground Level (AGL) by + subtracting the terrain height. + + This functions extracts the necessary variables from the NetCDF file + object in order to perform the calculation. + + Args: + + wrfin (:class:`netCDF4.Dataset`, :class:`Nio.NioFile`, or an \ + iterable): WRF-ARW NetCDF + data as a :class:`netCDF4.Dataset`, :class:`Nio.NioFile` + or an iterable sequence of the aforementioned types. + + timeidx (:obj:`int` or :data:`wrf.ALL_TIMES`, optional): The + desired time index. This value can be a positive integer, + negative integer, or + :data:`wrf.ALL_TIMES` (an alias for None) to return + all times in the file or sequence. The default is 0. + + method (:obj:`str`, optional): The aggregation method to use for + sequences. Must be either 'cat' or 'join'. + 'cat' combines the data along the Time dimension. + 'join' creates a new dimension for the file index. + The default is 'cat'. + + squeeze (:obj:`bool`, optional): Set to False to prevent dimensions + with a size of 1 from being automatically removed from the shape + of the output. Default is True. + + cache (:obj:`dict`, optional): A dictionary of (varname, ndarray) + that can be used to supply pre-extracted NetCDF variables to the + computational routines. It is primarily used for internal + purposes, but can also be used to improve performance by + eliminating the need to repeatedly extract the same variables + used in multiple diagnostics calculations, particularly when using + large sequences of files. + Default is None. + + meta (:obj:`bool`, optional): Set to False to disable metadata and + return :class:`numpy.ndarray` instead of + :class:`xarray.DataArray`. Default is True. + + _key (:obj:`int`, optional): A caching key. This is used for internal + purposes only. Default is None. + + units (:obj:`str`): The desired units. Refer to the :meth:`getvar` + product table for a list of available units for 'height_agl'. + Default is 'm'. + + Returns: + :class:`xarray.DataArray` or :class:`numpy.ndarray`: The + geopotential height. + If xarray is enabled and the *meta* parameter is True, then the result + will be a :class:`xarray.DataArray` object. Otherwise, the result will + be a :class:`numpy.ndarray` object with no metadata. + + """ + + return _get_geoht(wrfin, timeidx, method, squeeze, cache, meta, _key, + True, False) diff --git a/src/wrf/routines.py b/src/wrf/routines.py index 93727a5..3998209 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -8,7 +8,8 @@ from .g_cape import (get_2dcape, get_3dcape, get_cape2d_only, from .g_ctt import get_ctt from .g_dbz import get_dbz, get_max_dbz from .g_dewpoint import get_dp, get_dp_2m -from .g_geoht import get_geopt, get_height, get_stag_geopt, get_stag_height +from .g_geoht import (get_geopt, get_height, get_stag_geopt, get_stag_height, + get_height_agl) from .g_helicity import get_srh, get_uh from .g_latlon import get_lat, get_lon from .g_omega import get_omega @@ -78,6 +79,7 @@ _FUNC_MAP = {"cape2d": get_2dcape, "cloudfrac": get_cloudfrac, "geopt_stag": get_stag_geopt, "zstag": get_stag_height, + "height_agl": get_height_agl, # Diagnostics below are extracted from multi-product diagnostics "cape2d_only": get_cape2d_only, "cin2d_only": get_cin2d_only, @@ -143,6 +145,7 @@ _VALID_KARGS = {"cape2d": ["missing"], "mid_thresh", "high_thresh"], "geopt_stag": [], "zstag": ["msl", "units"], + "height_agl": ["units"], "cape2d_only": ["missing"], "cin2d_only": ["missing"], "lcl": ["missing"], diff --git a/test/ctt_test.py b/test/ctt_test.py index 6e97909..e5c625d 100644 --- a/test/ctt_test.py +++ b/test/ctt_test.py @@ -4,10 +4,12 @@ from matplotlib.cm import get_cmap import cartopy.crs as crs from cartopy.feature import NaturalEarthFeature -from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords +from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, + cartopy_ylim, latlon_coords) # Open the NetCDF file -ncfile = Dataset("/Users/ladwig/Documents/wrf_files/problem_files/cfrac_bug/wrfout_d02_1987-10-01_00:00:00") +ncfile = Dataset("/Users/ladwig/Documents/wrf_files/" + "problem_files/cfrac_bug/wrfout_d02_1987-10-01_00:00:00") # Get the sea level pressure ctt = getvar(ncfile, "ctt") @@ -19,20 +21,23 @@ lats, lons = latlon_coords(ctt) cart_proj = get_cartopy(ctt) # Create a figure -fig = plt.figure(figsize=(12,9)) +fig = plt.figure(figsize=(12, 9)) # Set the GeoAxes to the projection used by WRF ax = plt.axes(projection=cart_proj) # Download and add the states and coastlines -states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', +states = NaturalEarthFeature(category='cultural', scale='50m', + facecolor='none', name='admin_1_states_provinces_shp') ax.add_feature(states, linewidth=.5) ax.coastlines('50m', linewidth=0.8) -# Make the contour outlines and filled contours for the smoothed sea level pressure. +# Make the contour outlines and filled contours for the smoothed sea level +# pressure. plt.contour(to_np(lons), to_np(lats), to_np(ctt), 10, colors="black", transform=crs.PlateCarree()) -plt.contourf(to_np(lons), to_np(lats), to_np(ctt), 10, transform=crs.PlateCarree(), +plt.contourf(to_np(lons), to_np(lats), to_np(ctt), 10, + transform=crs.PlateCarree(), cmap=get_cmap("jet")) # Add a color bar diff --git a/test/generator_test.py b/test/generator_test.py index 6a45d89..8cc2816 100644 --- a/test/generator_test.py +++ b/test/generator_test.py @@ -1,17 +1,19 @@ -from __future__ import (absolute_import, division, print_function, unicode_literals) +from __future__ import (absolute_import, division, print_function, + unicode_literals) from wrf import getvar from netCDF4 import Dataset as nc -#ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00") + ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00") + def gen_seq(): wrfseq = [ncfile, ncfile, ncfile] for wrf in wrfseq: yield wrf - + + p_gen = getvar(gen_seq(), "P", method="join") print(p_gen) del p_gen - diff --git a/test/misc/extract_one_time.py b/test/misc/extract_one_time.py new file mode 100644 index 0000000..9645f98 --- /dev/null +++ b/test/misc/extract_one_time.py @@ -0,0 +1,40 @@ + +from netCDF4 import Dataset + +VARS_TO_KEEP = ("Times", "XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", + "XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", + "T2", "PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD", + "QGRAUP", "QRAIN", "QSNOW", "QICE", "MAPFAC_M", "MAPFAC_U", + "MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC", + "I_RAINC", "I_RAINNC") + +INPUT_FILE = "wrfout_d02_2005-08-28_00:00:00" +OUTPUT_FILE = "wrfout_problem_file" +DESIRED_TIME_INDEX = 0 + +with Dataset(INPUT_FILE) as infile, Dataset(OUTPUT_FILE, mode="w") as outfile: + # Copy the global attributes + outfile.setncatts(infile.__dict__) + + # Copy Dimensions + for name, dimension in infile.dimensions.items(): + if name != "Time": + dimsize = len(dimension) + outfile.createDimension(name, dimsize) + else: + outfile.createDimension(name, 1) + + # Copy Variables + for name, variable in infile.variables.iteritems(): + if name not in VARS_TO_KEEP: + continue + + outvar = outfile.createVariable(name, variable.datatype, + variable.dimensions, zlib=True) + + if len(variable.dimensions) > 1: + outvar[:] = variable[DESIRED_TIME_INDEX, :] + else: + outvar[:] = variable[DESIRED_TIME_INDEX] + + outvar.setncatts(variable.__dict__) diff --git a/test/misc/loop_and_fill_meta.py b/test/misc/loop_and_fill_meta.py new file mode 100644 index 0000000..7f46217 --- /dev/null +++ b/test/misc/loop_and_fill_meta.py @@ -0,0 +1,70 @@ +from __future__ import print_function, division + +import os +import numpy as np +from netCDF4 import Dataset +from wrf import getvar, ALL_TIMES, to_np +import xarray + +filename_list = ["/Users/ladwig/Documents/wrf_files/" + "wrf_vortex_single/wrfout_d02_2005-08-28_00:00:00", + "/Users/ladwig/Documents/wrf_files/wrf_vortex_single/" + "wrfout_d02_2005-08-28_03:00:00", + "/Users/ladwig/Documents/wrf_files/wrf_vortex_single/" + "wrfout_d02_2005-08-28_06:00:00", + "/Users/ladwig/Documents/wrf_files/wrf_vortex_single/" + "wrfout_d02_2005-08-28_09:00:00"] + +result_shape = (4, 1, 96, 96) + +# Let's get the first time so we can copy the metadata later +f = Dataset(filename_list[0]) +# By setting squeeze to False, you'll get all the dimension names. +z1 = getvar(f, "T2", 0, squeeze=False) +xlat = getvar(f, "XLAT", 0) +xlong = getvar(f, "XLONG", 0) + + +z_final = np.empty(result_shape, np.float32) + +# Modify this number if using more than 1 time per file +times_per_file = 1 + +data_times = [] +xtimes = [] +for timeidx in range(result_shape[0]): + # Compute the file index and the time index inside the file + fileidx = timeidx // times_per_file + file_timeidx = timeidx % times_per_file + + f = Dataset(filename_list[fileidx]) + z = getvar(f, "T2", file_timeidx) + t = getvar(f, "Times", file_timeidx) + xt = getvar(f, "xtimes", file_timeidx) + data_times.append(to_np(t)) + xtimes.append(to_np(xt)) + z_final[timeidx, :] = z[:] + f.close() + +# Let's make the metadata. Dimension names should copy easily if you set +# sqeeze to False, otherwise you can just set them yourself is a tuple of +# dimension names. Since you wanted +# to keep the bottom_top dimension for this 2D variable (which is normally +# removed), I'm doing this manually. +z_dims = ["Time", "bottom_top", "south_north", "west_east"] + +# Xarray doesn't copy coordinates easily (it always complains about shape +# mismatches), so do this manually +z_coords = {} +z_coords["Time"] = data_times +z_coords["XTIME"] = ("Time",), xtimes +z_coords["XLAT"] = ("south_north", "west_east"), xlat +z_coords["XLONG"] = ("south_north", "west_east"), xlong +z_name = "T2" + +# Attributes copy nicely +z_attrs = {} +z_attrs.update(z1.attrs) + +z_with_meta = xarray.DataArray(z_final, coords=z_coords, dims=z_dims, + attrs=z_attrs, name=z_name) diff --git a/test/misc/mocktest.py b/test/misc/mocktest.py index 211936f..7dd552a 100644 --- a/test/misc/mocktest.py +++ b/test/misc/mocktest.py @@ -6,39 +6,48 @@ try: except ImportError: from mock import Mock as MagicMock + class Mock(MagicMock): @classmethod def __getattr__(cls, name): return Mock() - -MOCK_MODULES = ["numpy", "numpy.ma", "xarray", "cartopy", + + +MOCK_MODULES = ["numpy", "numpy.ma", "xarray", "cartopy", "pandas", "matplotlib", "netCDF4", "mpl_toolkits.basemap", "wrf._wrffortran"] sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES) -consts = {"DEFAULT_FILL" : 9.9692099683868690E36, - "DEFAULT_FILL_INT8" : -127, - "DEFAULT_FILL_INT16" : -32767, - "DEFAULT_FILL_INT32" : -2147483647, - "DEFAULT_FILL_INT64" : -9223372036854775806, - "DEFAULT_FILL_FLOAT" : 9.9692099683868690E36, - "DEFAULT_FILL_DOUBLE" : 9.9692099683868690E36, - "fomp_sched_static" : 1, - "fomp_sched_dynamic" : 2, - "fomp_sched_guided" : 3, - "fomp_sched_auto" : 4} +consts = {"DEFAULT_FILL": 9.9692099683868690E36, + "DEFAULT_FILL_INT8": -127, + "DEFAULT_FILL_INT16": -32767, + "DEFAULT_FILL_INT32": -2147483647, + "DEFAULT_FILL_INT64": -9223372036854775806, + "DEFAULT_FILL_FLOAT": 9.9692099683868690E36, + "DEFAULT_FILL_DOUBLE": 9.9692099683868690E36, + "fomp_sched_static": 1, + "fomp_sched_dynamic": 2, + "fomp_sched_guided": 3, + "fomp_sched_auto": 4} + class MockWrfConstants(object): def __init__(self): self.__dict__ = consts - + + def mock_asscalar(val): - return float(val) + return float(val) + sys.modules["wrf._wrffortran"].wrf_constants = MockWrfConstants() sys.modules["wrf._wrffortran"].omp_constants = MockWrfConstants() sys.modules["numpy"].asscalar = mock_asscalar -import wrf -print (wrf.get_coord_pairs.__doc__) +try: + import wrf +except ImportError: + pass + +print(wrf.get_coord_pairs.__doc__) diff --git a/test/misc/projtest.py b/test/misc/projtest.py index 508b799..e1ac4d1 100644 --- a/test/misc/projtest.py +++ b/test/misc/projtest.py @@ -7,13 +7,15 @@ import matplotlib.cm as cm from netCDF4 import Dataset as NetCDF +from wrf import get_proj_params +from wrf.projection import getproj, RotatedLatLon, PolarStereographic PYNGL = True try: import Ngl except ImportError: PYNGL = False - + BASEMAP = True try: import mpl_toolkits.basemap @@ -25,19 +27,15 @@ try: from cartopy import crs, feature except ImportError: CARTOPY = False - -from wrf import get_proj_params -from wrf.projection import getproj, RotatedLatLon, PolarStereographic FILE_DIR = "/Users/ladwig/Documents/wrf_files/" -WRF_FILES = [ - join(FILE_DIR, "norway", "geo_em.d01.nc"), - join(FILE_DIR, "rotated_pole", "EAS_geo_em.d01.nc"), - join(FILE_DIR, "rotated_pole", "EUR_geo_em.d01.nc"), - join(FILE_DIR,"wrfout_d01_2016-02-25_18_00_00"), - join(FILE_DIR, "wrfout_d01_2008-09-29_23-30-00"), - join(FILE_DIR, "wrfout_d01_2010-06-13_21:00:00")] +WRF_FILES = [join(FILE_DIR, "norway", "geo_em.d01.nc"), + join(FILE_DIR, "rotated_pole", "EAS_geo_em.d01.nc"), + join(FILE_DIR, "rotated_pole", "EUR_geo_em.d01.nc"), + join(FILE_DIR, "wrfout_d01_2016-02-25_18_00_00"), + join(FILE_DIR, "wrfout_d01_2008-09-29_23-30-00"), + join(FILE_DIR, "wrfout_d01_2010-06-13_21:00:00")] def nz_proj(): @@ -45,65 +43,68 @@ def nz_proj(): [-32.669853, -32.669853]]) lons = np.array([[163.839595, -179.693502], [163.839595, -179.693502]]) - - params = {"MAP_PROJ" : 6, - "CEN_LAT" : -41.814869, - "CEN_LON" : 179.693502, - "TRUELAT1" : 0, + + params = {"MAP_PROJ": 6, + "CEN_LAT": -41.814869, + "CEN_LON": 179.693502, + "TRUELAT1": 0, "TRUELAT2": 0, - "MOAD_CEN_LAT" : -41.814869, - "STAND_LON" : 180.0 - 179.693502, - "POLE_LAT" : 48.185131, - "POLE_LON" : 0.0} - + "MOAD_CEN_LAT": -41.814869, + "STAND_LON": 180.0 - 179.693502, + "POLE_LAT": 48.185131, + "POLE_LON": 0.0} + return lats, lons, RotatedLatLon(lats=lats, lons=lons, **params) + def argentina_proj(): lats = np.array([[-57.144064, -57.144064], [-21.154470, -21.154470]]) lons = np.array([[-86.893797, -37.089724], [-86.893797, -37.089724]]) - - params = {"MAP_PROJ" : 6, - "CEN_LAT" : -39.222954, - "CEN_LON" : -65.980109, - "TRUELAT1" : 0, + + params = {"MAP_PROJ": 6, + "CEN_LAT": -39.222954, + "CEN_LON": -65.980109, + "TRUELAT1": 0, "TRUELAT2": 0, - "MOAD_CEN_LAT" : -39.222954, - "STAND_LON" : 180.0 - -65.980109, - "POLE_LAT" : 90 + -39.222954, - "POLE_LON" : 0.0} - + "MOAD_CEN_LAT": -39.222954, + "STAND_LON": 180.0 - -65.980109, + "POLE_LAT": 90 + -39.222954, + "POLE_LON": 0.0} + return lats, lons, RotatedLatLon(lats=lats, lons=lons, **params) + def south_polar_proj(): - lats = np.array([[-30.0,-30.0], - [-30.0,-30.0]]) + lats = np.array([[-30.0, -30.0], + [-30.0, -30.0]]) lons = np.array([[-120, 60], [-120, 60]]) - - params = {"MAP_PROJ" : 2, - "CEN_LAT" : -90.0, - "CEN_LON" : 0, - "TRUELAT1" : -10.0, - "MOAD_CEN_LAT" : -90.0, - "STAND_LON" : 0} - + + params = {"MAP_PROJ": 2, + "CEN_LAT": -90.0, + "CEN_LON": 0, + "TRUELAT1": -10.0, + "MOAD_CEN_LAT": -90.0, + "STAND_LON": 0} + return lats, lons, PolarStereographic(lats=lats, lons=lons, **params) + def north_polar_proj(): - lats = np.array([[30.0,30.0], - [30.0,30.0]]) + lats = np.array([[30.0, 30.0], + [30.0, 30.0]]) lons = np.array([[-45, 140], [-45, 140]]) - - params = {"MAP_PROJ" : 2, - "CEN_LAT" : 90.0, - "CEN_LON" : 10, - "TRUELAT1" : 10.0, - "MOAD_CEN_LAT" : 90.0, - "STAND_LON" : 10} - + + params = {"MAP_PROJ": 2, + "CEN_LAT": 90.0, + "CEN_LON": 10, + "TRUELAT1": 10.0, + "MOAD_CEN_LAT": 90.0, + "STAND_LON": 10} + return lats, lons, PolarStereographic(lats=lats, lons=lons, **params) @@ -112,21 +113,23 @@ def dateline_rot_proj(): [71.717521, 71.717521]]) lons = np.array([[170.332771, -153.456292], [170.332771, -153.456292]]) - - params = {"MAP_PROJ" : 6, - "CEN_LAT" : 66.335764, - "CEN_LON" : -173.143792, - "TRUELAT1" : 0, + + params = {"MAP_PROJ": 6, + "CEN_LAT": 66.335764, + "CEN_LON": -173.143792, + "TRUELAT1": 0, "TRUELAT2": 0, - "MOAD_CEN_LAT" : 66.335764, - "STAND_LON" : 173.143792, - "POLE_LAT" : 90.0 - 66.335764, - "POLE_LON" : 180.0} + "MOAD_CEN_LAT": 66.335764, + "STAND_LON": 173.143792, + "POLE_LAT": 90.0 - 66.335764, + "POLE_LON": 180.0} return lats, lons, RotatedLatLon(lats=lats, lons=lons, **params) + class WRFProjTest(ut.TestCase): longMessage = True + def make_test(wrf_file=None, fixed_case=None): if wrf_file is not None: ncfile = NetCDF(wrf_file) @@ -144,75 +147,76 @@ def make_test(wrf_file=None, fixed_case=None): elif fixed_case == "north_polar": lats, lons, proj = north_polar_proj() elif fixed_case == "dateline_rot": - lats,lons,proj = dateline_rot_proj() - - print ("wrf proj4: {}".format(proj.proj4())) + lats, lons, proj = dateline_rot_proj() + + print("wrf proj4: {}".format(proj.proj4())) if PYNGL: # PyNGL plotting wks_type = bytes("png") - wks = Ngl.open_wks(wks_type,bytes("pyngl_{}".format(name_suffix))) + wks = Ngl.open_wks(wks_type, bytes("pyngl_{}".format(name_suffix))) mpres = proj.pyngl() - map = Ngl.map(wks,mpres) - + map = Ngl.map(wks, mpres) + Ngl.delete_wks(wks) - - if BASEMAP: + + if BASEMAP: # Basemap plotting - fig = plt.figure(figsize=(10,10)) - ax = fig.add_axes([0.1,0.1,0.8,0.8]) - + fig = plt.figure(figsize=(10, 10)) + ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) + # Define and plot the meridians and parallels min_lat = np.amin(lats) max_lat = np.amax(lats) min_lon = np.amin(lons) max_lon = np.amax(lons) - - parallels = np.arange(np.floor(min_lat), np.ceil(max_lat), + + parallels = np.arange(np.floor(min_lat), np.ceil(max_lat), (max_lat - min_lat)/5.0) - meridians = np.arange(np.floor(min_lon), np.ceil(max_lon), + meridians = np.arange(np.floor(min_lon), np.ceil(max_lon), (max_lon - min_lon)/5.0) - + bm = proj.basemap() bm.drawcoastlines(linewidth=.5) - #bm.drawparallels(parallels,labels=[1,1,1,1],fontsize=10) - #bm.drawmeridians(meridians,labels=[1,1,1,1],fontsize=10) - print ("basemap proj4: {}".format(bm.proj4string)) + # bm.drawparallels(parallels,labels=[1,1,1,1],fontsize=10) + # bm.drawmeridians(meridians,labels=[1,1,1,1],fontsize=10) + print("basemap proj4: {}".format(bm.proj4string)) plt.savefig("basemap_{}.png".format(name_suffix)) plt.close(fig) - + if CARTOPY: # Cartopy plotting - fig = plt.figure(figsize=(10,10)) - ax = plt.axes([0.1,0.1,0.8,0.8], projection=proj.cartopy()) - print ("cartopy proj4: {}".format(proj.cartopy().proj4_params)) - + fig = plt.figure(figsize=(10, 10)) + ax = plt.axes([0.1, 0.1, 0.8, 0.8], projection=proj.cartopy()) + print("cartopy proj4: {}".format(proj.cartopy().proj4_params)) + ax.coastlines('50m', linewidth=0.8) - #print proj.x_extents() - #print proj.y_extents() + # print proj.x_extents() + # print proj.y_extents() ax.set_xlim(proj.cartopy_xlim()) ax.set_ylim(proj.cartopy_ylim()) ax.gridlines() plt.savefig("cartopy_{}.png".format(name_suffix)) plt.close(fig) + if __name__ == "__main__": for wrf_file in WRF_FILES: test_func = make_test(wrf_file=wrf_file) setattr(WRFProjTest, "test_proj", test_func) - + test_func2 = make_test(fixed_case="south_rot") setattr(WRFProjTest, "test_south_rot", test_func2) - + test_func3 = make_test(fixed_case="arg_rot") setattr(WRFProjTest, "test_arg_rot", test_func3) - + test_func4 = make_test(fixed_case="south_polar") setattr(WRFProjTest, "test_south_polar", test_func4) - + test_func5 = make_test(fixed_case="north_polar") setattr(WRFProjTest, "test_north_polar", test_func5) - + test_func6 = make_test(fixed_case="dateline_rot") setattr(WRFProjTest, "test_dateline_rot", test_func6) - - ut.main() \ No newline at end of file + + ut.main() diff --git a/test/misc/quiver_test.py b/test/misc/quiver_test.py new file mode 100644 index 0000000..557167e --- /dev/null +++ b/test/misc/quiver_test.py @@ -0,0 +1,58 @@ +from netCDF4 import Dataset +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.cm import get_cmap +import cartopy.crs as crs +from cartopy.feature import NaturalEarthFeature + +from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy, + cartopy_xlim, cartopy_ylim) + +# Open the NetCDF file +ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + +# Extract the pressure and wind variables + +p = getvar(ncfile, "pressure") +# Note: This is m/s. +ua = getvar(ncfile, "ua") +va = getvar(ncfile, "va") + +# Interpolate u, and v winds to 950 hPa +u_950 = interplevel(ua, p, 950) +v_950 = interplevel(va, p, 950) + +# Get the lat/lon coordinates +lats, lons = latlon_coords(u_950) + +# Get the map projection information +cart_proj = get_cartopy(u_950) + +# Create the figure +fig = plt.figure(figsize=(12, 9)) +ax = plt.axes(projection=cart_proj) + +# Download and add the states and coastlines +states = NaturalEarthFeature(category='cultural', scale='50m', + facecolor='none', + name='admin_1_states_provinces_shp') +ax.add_feature(states, linewidth=0.5) +ax.coastlines('50m', linewidth=0.8) + +# Add the 950 hPa wind barbs, only plotting every 'thin'ed barb. Adjust thin +# as needed. Also, no scaling is done for the arrows, so you might need to +# mess with the scale argument. +thin = 50 +plt.quiver(to_np(lons[::thin, ::thin]), to_np(lats[::thin, ::thin]), + to_np(u_950[::thin, ::thin]), to_np(v_950[::thin, ::thin]), + transform=crs.PlateCarree()) + +# Set the map bounds +ax.set_xlim(cartopy_xlim(u_950)) +ax.set_ylim(cartopy_ylim(v_950)) + +ax.gridlines() + +plt.title("Arrows") + +plt.show() diff --git a/test/misc/reduce_file.py b/test/misc/reduce_file.py index 35e61d9..d7e0b08 100644 --- a/test/misc/reduce_file.py +++ b/test/misc/reduce_file.py @@ -13,22 +13,22 @@ for att_name in in_nc.ncattrs(): # Copy Dimensions, but modify the vertical dimensions for dimname, dim in in_nc.dimensions.iteritems(): out_nc.createDimension(dimname, len(dim)) - + # Copy Variables and their Attributes, using the reduced vertical dimension for varname, var in in_nc.variables.iteritems(): if varname in ("T2", "XLAT", "XLONG", "XTIME"): datatype = var.datatype dimensions = var.dimensions shape = var.shape - + new_shape = shape - + new_var = out_nc.createVariable(varname, datatype, dimensions) - + new_var[:] = var[:] - + for att_name in var.ncattrs(): setattr(new_var, att_name, getattr(var, att_name)) - + in_nc.close() -out_nc.close() \ No newline at end of file +out_nc.close() diff --git a/test/misc/snippet.py b/test/misc/snippet.py index 014810f..d3df43a 100644 --- a/test/misc/snippet.py +++ b/test/misc/snippet.py @@ -1,30 +1,29 @@ - import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap + def main(): - bm = Basemap(projection = "rotpole", - o_lat_p = 36.0, - o_lon_p = 180.0, - llcrnrlat = -10.590603, - urcrnrlat = 46.591976, - llcrnrlon = -139.08585, - urcrnrlon = 22.661009, - lon_0 = -106.0, - rsphere = 6370000, - resolution = 'l') - - fig = plt.figure(figsize=(8,8)) - ax = fig.add_axes([0.1,0.1,0.8,0.8]) - + bm = Basemap(projection="rotpole", + o_lat_p=36.0, + o_lon_p=180.0, + llcrnrlat=-10.590603, + urcrnrlat=46.591976, + llcrnrlon=-139.08585, + urcrnrlon=22.661009, + lon_0=-106.0, + rsphere=6370000, + resolution='l') + + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) + bm.drawcoastlines(linewidth=.5) - - print bm.proj4string - + + print(bm.proj4string) + plt.savefig("basemap_map.png") plt.close(fig) - - + if __name__ == "__main__": main() diff --git a/test/varcache.py b/test/misc/varcache.py similarity index 68% rename from test/varcache.py rename to test/misc/varcache.py index f8ba15d..54ca5ed 100644 --- a/test/varcache.py +++ b/test/misc/varcache.py @@ -4,33 +4,37 @@ import time from netCDF4 import Dataset from wrf import getvar, ALL_TIMES, extract_vars -wrf_filenames = ["/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-28_00:00:00", - "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-28_12:00:00", - "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-29_00:00:00"] +wrf_filenames = ["/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/" + "moving_nest/wrfout_d02_2005-08-28_00:00:00", + "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/" + "moving_nest/wrfout_d02_2005-08-28_12:00:00", + "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/" + "moving_nest/wrfout_d02_2005-08-29_00:00:00"] wrfin = [Dataset(x) for x in wrf_filenames] -my_cache = extract_vars(wrfin, ALL_TIMES, ("P", "PB", "PH", "PHB", "T", "QVAPOR", "HGT", "U", "V", "W", "PSFC")) +my_cache = extract_vars(wrfin, ALL_TIMES, + ("P", "PB", "PH", "PHB", "T", "QVAPOR", "HGT", "U", + "V", "W", "PSFC")) start = time.time() -for var in ("avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", - "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", - "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", - "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", +for var in ("avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", + "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", + "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", + "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"): v = getvar(wrfin, var, ALL_TIMES) end = time.time() -print ("Time taken without variable cache: ", (end-start)) +print("Time taken without variable cache: ", (end-start)) start = time.time() -for var in ("avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", - "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", - "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", - "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", +for var in ("avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", + "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", + "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", + "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"): v = getvar(wrfin, var, ALL_TIMES, cache=my_cache) end = time.time() -print ("Time taken with variable cache: ", (end-start)) - +print("Time taken with variable cache: ", (end-start)) diff --git a/test/viewtest.py b/test/misc/viewtest.py similarity index 51% rename from test/viewtest.py rename to test/misc/viewtest.py index 5dab455..f9a7064 100644 --- a/test/viewtest.py +++ b/test/misc/viewtest.py @@ -1,26 +1,25 @@ +# Test snippet for f2py import numpy as np import wrf._wrffortran errlen = int(wrf._wrffortran.constants.errlen) - -a = np.ones((3,3,3)) -b = np.zeros((3,3,3,3)) +a = np.ones((3, 3, 3)) +b = np.zeros((3, 3, 3, 3)) c = np.zeros(errlen, "c") errstat = np.array(0) errmsg = np.zeros(errlen, "c") c[:] = "Test String" - for i in xrange(2): - outview = b[i,:] + outview = b[i, :] outview = outview.T - q = wrf._wrffortran.testfunc(a,outview,c,errstat=errstat,errmsg=errmsg) - print errstat - + q = wrf._wrffortran.testfunc(a, outview, c, errstat=errstat, + errmsg=errmsg) + print(errstat) str_bytes = (bytes(char).decode("utf-8") for char in errmsg[:]) -print repr(errmsg) -print "".join(str_bytes).strip() \ No newline at end of file +print(repr(errmsg)) +print("".join(str_bytes).strip()) diff --git a/test/misc/wps.py b/test/misc/wps.py new file mode 100644 index 0000000..b2a80ea --- /dev/null +++ b/test/misc/wps.py @@ -0,0 +1,230 @@ +# Hastily made script to read WPS intermediate files +from __future__ import print_function +import struct + +import numpy as np + +# Number of bytes used at the start and end of a fortran record to +# indicate the record size +SIZE_BYTES = 4 + + +class WPSData(object): + def __init__(self, ifv=None, hdate=None, xfcst=None, map_source=None, + field=None, units=None, desc=None, xlvl=None, nx=None, + ny=None, iproj=None, startloc=None, startlat=None, + startlon=None, deltalat=None, deltalon=None, nlats=None, + dx=None, dy=None, xlonc=None, truelat1=None, truelat2=None, + earth_radius=None, is_wind_earth_rel=None, slab=None): + + self.ifv = ifv + self.hdate = hdate + self.xfcst = xfcst + self.map_source = map_source + self.field = field + self.units = units + self.desc = desc + self.xlvl = xlvl + self.nx = nx + self.ny = ny + self.iproj = iproj + self.startloc = startloc + self.startlat = startlat + self.startlon = startlon + self.deltalat = deltalat + self.deltalon = deltalon + self.nlats = nlats + self.dx = dx + self.dy = dy + self.xlonc = xlonc + self.truelat1 = truelat1 + self.truelat2 = truelat2 + self.earth_radius = earth_radius + self.is_wind_earth_rel = is_wind_earth_rel + self.slab = slab + + +def _parse_record1(data): + result = {} + result["ifv"] = struct.unpack(">i", data) + + return result + + +def _parse_record2(data): + result = {} + parsed = struct.unpack(">24sf32s9s25s46sfiii", data) + result["hdate"] = parsed[0] + result["xfcst"] = parsed[1] + result["map_source"] = parsed[2] + result["field"] = parsed[3] + result["units"] = parsed[4] + result["desc"] = parsed[5] + result["xlvl"] = parsed[6] + result["nx"] = parsed[7] + result["ny"] = parsed[8] + result["iproj"] = parsed[9] + + return result + + +def _parse_record3(data, iproj): + result = {} + if iproj == 0: + fmt = ">8sfffff" + parsed = struct.unpack(fmt, data) + result["startloc"] = parsed[0] + result["startlat"] = parsed[1] + result["startlon"] = parsed[2] + result["deltalat"] = parsed[3] + result["deltalon"] = parsed[4] + result["earth_radius"] = parsed[5] + elif iproj == 1: + fmt = ">8sffffff" + parsed = struct.unpack(fmt, data) + result["startloc"] = parsed[0] + result["startlat"] = parsed[1] + result["startlon"] = parsed[2] + result["dx"] = parsed[3] + result["dy"] = parsed[4] + result["truelat1"] = parsed[5] + result["earth_radius"] = parsed[6] + elif iproj == 3: + fmt = ">8sffffffff" + parsed = struct.unpack(fmt, data) + result["startloc"] = parsed[0] + result["startlat"] = parsed[1] + result["startlon"] = parsed[2] + result["dx"] = parsed[3] + result["dy"] = parsed[4] + result["xlonc"] = parsed[5] + result["truelat1"] = parsed[6] + result["truelat2"] = parsed[7] + result["earth_radius"] = parsed[8] + elif iproj == 4: + fmt = ">8sfffff" + parsed = struct.unpack(fmt, data) + result["startloc"] = parsed[0] + result["startlat"] = parsed[1] + result["startlon"] = parsed[2] + result["nlats"] = parsed[3] + result["deltalon"] = parsed[4] + result["earth_radius"] = parsed[5] + elif iproj == 5: + fmt = ">8sfffffff" + parsed = struct.unpack(fmt, data) + result["startloc"] = parsed[0] + result["startlat"] = parsed[1] + result["startlon"] = parsed[2] + result["dx"] = parsed[3] + result["dy"] = parsed[4] + result["xlonc"] = parsed[5] + result["truelat1"] = parsed[6] + result["earth_radius"] = parsed[7] + + return result + + +def _parse_record4(data): + result = {} + result["is_wind_earth_rel"] = struct.unpack(">i", data) + + return result + + +def _parse_record5(data, nx, ny): + result = {} + + size = nx * ny + fmt = ">{}f".format(size) + parsed = struct.unpack(fmt, data) + + arr = np.array(parsed, dtype=np.float32) + result["slab"] = arr.reshape((nx, ny), order="F") + + return result + + +_PARSE_MAP = {0: _parse_record1, + 1: _parse_record2, + 2: _parse_record3, + 3: _parse_record4, + 4: _parse_record5} + + +def parse_record(record_idx, data, iproj=None, nx=None, ny=None): + + if record_idx == 0 or record_idx == 1 or record_idx == 3: + return _PARSE_MAP[record_idx](data) + elif record_idx == 2: + return _PARSE_MAP[record_idx](data, iproj) + elif record_idx == 4: + return _PARSE_MAP[record_idx](data, nx, ny) + + +def read_wps(wps_file, field=None): + with open(wps_file, "rb") as f: + wps_params = {} + record_list = [] + raw_data = f.read() + + record_size_idx = 0 + end_of_file_idx = len(raw_data) - 1 + + while True: + iproj = None + nx = None + ny = None + keep_record = True + for record_idx in range(5): + # Each record has the size (in SIZE_BYTES bytes) at the + # start and end of each record. This might be compiler + # dependent though, so this might need to be modified. + # Also, the WPS files are stored big endian. + + record_size = struct.unpack( + ">i", + raw_data[record_size_idx: record_size_idx + SIZE_BYTES]) + record_start = record_size_idx + SIZE_BYTES + record_end = record_start + record_size[0] + record_data = raw_data[record_start:record_end] + + parsed_record = parse_record(record_idx, record_data, iproj, + nx, ny) + + try: + field_name = parsed_record["field"].strip() + except KeyError: + pass + else: + if field is not None: + if field_name != field: + keep_record = False + + try: + iproj = parsed_record["iproj"] + except KeyError: + pass + + try: + nx = parsed_record["nx"] + except KeyError: + pass + + try: + ny = parsed_record["ny"] + except KeyError: + pass + + wps_params.update(parsed_record) + + record_size_idx = record_end + SIZE_BYTES + + if keep_record: + record_list.append(WPSData(**wps_params)) + + # Repeat for all record slabs + if record_end + SIZE_BYTES > end_of_file_idx: + break + + return record_list diff --git a/test/listBug.ncl b/test/ncl/listBug.ncl similarity index 100% rename from test/listBug.ncl rename to test/ncl/listBug.ncl diff --git a/test/ncl_get_var.ncl b/test/ncl/ncl_get_var.ncl similarity index 100% rename from test/ncl_get_var.ncl rename to test/ncl/ncl_get_var.ncl diff --git a/test/ncl/ncl_vertcross.ncl b/test/ncl/ncl_vertcross.ncl new file mode 100644 index 0000000..5b9991a --- /dev/null +++ b/test/ncl/ncl_vertcross.ncl @@ -0,0 +1,92 @@ +input_file = addfile("/Users/ladwig/Documents/wrf_files/wrfout_d02_2010-06-13_21:00:00.nc", "r") + +z = wrf_user_getvar(input_file, "z", 0) ; grid point height +p = wrf_user_getvar(input_file, "pressure", 0) ; total pressure + +dimsz = dimsizes(z) +pivot = (/ (dimsz(2)-1)/2, (dimsz(1)-1)/2 /) ; pivot point is center of domain + +; For the new cross section routine +xopt = True +xopt@use_pivot = True +xopt@angle = 45.0 +;xopt@levels = +;xopt@latlon = +xopt@file_handle = input_file +;xopt@timeidx = +xopt@linecoords = True + +ht_vertcross = wrf_user_vertcross(z, p, pivot, xopt) + +printVarSummary(ht_vertcross) +print(min(ht_vertcross@lats)) +print(min(ht_vertcross@lons)) +print(max(ht_vertcross@lats)) +print(max(ht_vertcross@lons)) + + +xopt@use_pivot = False +xopt@angle = 0.0 +;xopt@levels = +xopt@latlon = True +xopt@file_handle = input_file +xopt@timeidx = 0 +xopt@linecoords = True + +loc_param = (/-104.3632, 32.8562, -95.15308, 40.06575 /) ; pivot point is center of domain +ht_vertcross2 = wrf_user_vertcross(z, p, loc_param, xopt) + +printVarSummary(ht_vertcross2) +print(min(ht_vertcross2@lats)) +print(min(ht_vertcross2@lons)) +print(max(ht_vertcross2@lats)) +print(max(ht_vertcross2@lons)) + +print(ht_vertcross2@lats(190)) +print(ht_vertcross2@lons(190)) + +xopt@use_pivot = True +xopt@angle = 45.0 +;xopt@levels = +xopt@latlon = True +xopt@file_handle = input_file +xopt@timeidx = 0 +xopt@linecoords = True + +loc_param := (/-99.98572, 36.54949 /) ; pivot point is center of domain +ht_vertcross3 = wrf_user_vertcross(z, p, loc_param, xopt) + +printVarSummary(ht_vertcross3) +print(min(ht_vertcross3@lats)) +print(min(ht_vertcross3@lons)) +print(max(ht_vertcross3@lats)) +print(max(ht_vertcross3@lons)) + + +xopt@use_pivot = True +xopt@angle = 45.0 +xopt@levels = (/1000., 850., 700., 500., 250. /) +xopt@latlon = True +xopt@file_handle = input_file +xopt@timeidx = 0 +xopt@linecoords = True + +loc_param := (/-99.98572, 36.54949 /) ; pivot point is center of domain +ht_vertcross4 = wrf_user_vertcross(z, p, loc_param, xopt) + +printVarSummary(ht_vertcross4) +print(min(ht_vertcross4@lats)) +print(min(ht_vertcross4@lons)) +print(max(ht_vertcross4@lats)) +print(max(ht_vertcross4@lons)) + +o = True +o@returnInt = False +o@useTime = 0 +l = wrf_user_ll_to_xy(input_file, -99.98572, 36.54949, o) +print(l) + + +l1 = wrf_user_xy_to_ll(input_file, l(1), l(0), o) +print(l1) + diff --git a/test/ncl/refl10_cross.ncl b/test/ncl/refl10_cross.ncl new file mode 100644 index 0000000..c90861e --- /dev/null +++ b/test/ncl/refl10_cross.ncl @@ -0,0 +1,81 @@ + a = addfile("wrfout_d03_2017-04-03_06:00:00_ctrl","r") + + time = 0 + refl_10cm = wrf_user_getvar(a,"REFL_10CM",time) + z = wrf_user_getvar(a, "z", time) + lat = wrf_user_getvar(a, "lat", time) + lon = wrf_user_getvar(a, "lon", time) + + ; convert the lat/lon to x,y + start_lat = 20.9 + start_lon = 92.5 + end_lat = 29.2 + end_lon = 92.5 + + opt = True + + start_ij = wrf_user_ll_to_ij(a, start_lon, start_lat, opt) + start_ij = start_ij - 1 + + end_ij = wrf_user_ll_to_ij(a, end_lon, end_lat, opt) + end_ij = end_ij - 1 + + start_end = (/start_ij(0), start_ij(1), end_ij(0), end_ij(1)/) + + lat_line = wrf_user_intrp2d(lat,start_end,0.0,True) + nlat = dimsizes(lat_line) + + lon_line = wrf_user_intrp2d(lon,start_end,0.0,True) + + refl_cross = wrf_user_intrp3d(refl_10cm,z,"v",start_end,0.,True) + + ; Need to make a vertical coordinate by using the same code as the + ; cross section + + ; Currently, the vertical coordinate is not set, so let's do it + ; manually. This will be fixed in the next version of NCL. + ; If you want to set these levels yourself, you'll need to copy the + ; code I sent before and manually set the levels in the cross section + ; routine, then do it again here. + + z_max = max(z) + z_min = 0. + dz = 0.01 * z_max + nlevels = tointeger( z_max/dz ) + z_var2d = new( (/nlevels/), typeof(z)) + z_var2d(0) = z_min + + do i=1, nlevels-1 + z_var2d(i) = z_var2d(0)+i*dz + end do + + refl_cross&Vertical = z_var2d + + wks = gsn_open_wks("png","cross") + cmap := read_colormap_file("BlAqGrYeOrReVi200") + cmap(0,:) = (/0,0,0,0/) ; make first color fully transparent + + resx = True + resx@gsnMaximize = True + resx@lbLabelAutoStride = True ; default v6.1.0 + + resx@cnFillOn = True ; turn on color fill + resx@cnLinesOn = False ; turn lines on/off ; True is default + resx@cnLineLabelsOn = False ; turn line labels on/off ; True is default + resx@cnFillPalette = cmap + nLabels = 8 ; arbitrary + resx@tmXBLabels = new(nLabels,"string") + resx@tmXBMode = "Explicit" + + resx@tmXBValues := toint(fspan(0,nlat-1,nLabels)) + do i=0,nLabels-1 + x = lon_line(i) + y = lat_line(i) + resx@tmXBLabels(i) = sprintf("%5.1f", y)+"~C~"+sprintf("%5.1f", x) + end do + + resx@tiMainString = "Full South-North Grid Line X-Section" + + + plot1 = gsn_csm_contour(wks, refl_cross, resx ) + \ No newline at end of file diff --git a/test/ncl/rotated_test.ncl b/test/ncl/rotated_test.ncl new file mode 100644 index 0000000..9eb6f94 --- /dev/null +++ b/test/ncl/rotated_test.ncl @@ -0,0 +1,26 @@ + +;---Open file and calculate slp. +a = addfile("/Users/ladwig/Documents/wrf_files/rotated_pole_test.nc","r") + +t2 = wrf_user_getvar(a,"T2",0) + +;---Start the graphics +wks = gsn_open_wks("x11","wrf") + +;---Set some graphical resources +res = True +res@gsnMaximize = True +res@cnFillOn = True +res@tfDoNDCOverlay = True ; This is necessary if you don't + ; set sfXArray/sfYArray + +;---Add additional map resources +res = wrf_map_resources(a,res) + +;---Change some of the resources that were set (these were set to "gray") +res@mpGeophysicalLineColor = "black" +res@mpNationalLineColor = "black" +res@mpUSStateLineColor = "black" + +plot = gsn_csm_contour_map(wks,t2,res) + diff --git a/test/ncl/test_this.ncl b/test/ncl/test_this.ncl new file mode 100644 index 0000000..8421690 --- /dev/null +++ b/test/ncl/test_this.ncl @@ -0,0 +1,21 @@ +ifils = systemfunc ("ls /Users/ladwig/Documents/wrf_files/wrf_vortex_single/wrfout_d02_*") +print(ifils) +a = addfiles(ifils, "r") +;a = addfile("/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc", "r") + +lats := (/22.0, 25.0, 27.0 /) +lons := (/-90.0, -87.5, -83.75 /) + +opt = True +opt@useTime = -1 +opt@returnI = False +xy = wrf_user_ll_to_xy(a, lons, lats, opt) + +print(xy) + +x_s = (/10, 50, 90 /) +y_s = (/10, 50, 90 /) + +ll = wrf_user_xy_to_ll(a, x_s, y_s, opt) +print(ll) + diff --git a/test/test_vinterp.ncl b/test/ncl/test_vinterp.ncl similarity index 100% rename from test/test_vinterp.ncl rename to test/ncl/test_vinterp.ncl diff --git a/test/ncl/wrf_user_vertcross.ncl b/test/ncl/wrf_user_vertcross.ncl new file mode 100644 index 0000000..54df159 --- /dev/null +++ b/test/ncl/wrf_user_vertcross.ncl @@ -0,0 +1,416 @@ +;-------------------------------------------------------------------------------- + +undef("wrf_user_ll_to_xy") +function wrf_user_ll_to_xy( file_handle, longitude:numeric, latitude:numeric, \ + opts_args:logical ) + +; This is the same as wrf_user_ll_to_ij, but returns 0-based indexes + +begin +; +; As of NCL V6.0.0, wrf_user_ll_to_ij can now handle a file +; or a list of files. +; + if(typeof(file_handle).eq."file") then + ISFILE = True + nc_file = file_handle + elseif(typeof(file_handle).eq."list") then + ISFILE = False + nc_file = file_handle[0] + else + print("wrf_user_ll_to_xy: Error: the first argument must be a file or a list of files opened with addfile or addfiles") + return + end if + + opts = opts_args + useT = get_res_value(opts,"useTime",0) + returnI= get_res_value(opts,"returnInt",True) + + res = True + res@MAP_PROJ = nc_file@MAP_PROJ + res@TRUELAT1 = nc_file@TRUELAT1 + res@TRUELAT2 = nc_file@TRUELAT2 + res@STAND_LON = nc_file@STAND_LON + res@DX = nc_file@DX + res@DY = nc_file@DY + + if (res@MAP_PROJ .eq. 6) then + res@POLE_LAT = nc_file@POLE_LAT + res@POLE_LON = nc_file@POLE_LON + res@LATINC = (res@DY*360.)/2.0/3.141592653589793/6370000. + res@LONINC = (res@DX*360.)/2.0/3.141592653589793/6370000. + else + res@POLE_LAT = 90.0 + res@POLE_LON = 0.0 + res@LATINC = 0.0 + res@LONINC = 0.0 + end if + + if(isfilevar(nc_file,"XLAT")) + if(ISFILE) then + XLAT = nc_file->XLAT(useT,:,:) + XLONG = nc_file->XLONG(useT,:,:) + else + XLAT = file_handle[:]->XLAT + XLONG = file_handle[:]->XLONG + end if + else + if(ISFILE) then + XLAT = nc_file->XLAT_M(useT,:,:) + XLONG = nc_file->XLONG_M(useT,:,:) + else + XLAT = file_handle[:]->XLAT_M + XLONG = file_handle[:]->XLONG_M + end if + end if + + + if(dimsizes(dimsizes(XLAT)).eq.2) then +; Rank 2 + res@REF_LAT = XLAT(0,0) + res@REF_LON = XLONG(0,0) + else +; Rank 3 + res@REF_LAT = XLAT(useT,0,0) + res@REF_LON = XLONG(useT,0,0) + end if + res@KNOWNI = 1.0 + res@KNOWNJ = 1.0 + + loc = wrf_ll_to_ij (longitude, latitude, res) + loc = loc - 1 + + if (dimsizes(dimsizes(loc)) .eq. 1) then + loc!0 = "x_y" + elseif (dimsizes(dimsizes(loc)) .eq. 2) then + loc!0 = "x_y" + loc!1 = "idx" + else ; Not currently supported + loc!0 = "x_y" + loc!1 = "domain_idx" + loc!2 = "idx" + end if + + if ( returnI ) then + loci = new(dimsizes(loc),integer) + ;loci@_FillValue = default_fillvalue("integer") ; was -999 + loci = tointeger(loc + .5) + loci!0 = loc!0 + return(loci) + else + return(loc) + end if + + +end + +;-------------------------------------------------------------------------------- + +undef("wrf_user_xy_to_ll") +function wrf_user_xy_to_ll( file_handle, x:numeric, y:numeric, \ + opts_args:logical ) + +begin +; +; As of NCL V6.0.0, wrf_user_ll_to_ij can now handle a file +; or a list of files. +; + if(typeof(file_handle).eq."file") then + ISFILE = True + nc_file = file_handle + elseif(typeof(file_handle).eq."list") then + ISFILE = False + nc_file = file_handle[0] + else + print("wrf_user_xy_to_ll: Error: the first argument must be a file or a list of files opened with addfile or addfiles") + return + end if + + opts = opts_args + useT = get_res_value(opts,"useTime",0) + + res = True + res@MAP_PROJ = nc_file@MAP_PROJ + res@TRUELAT1 = nc_file@TRUELAT1 + res@TRUELAT2 = nc_file@TRUELAT2 + res@STAND_LON = nc_file@STAND_LON + res@DX = nc_file@DX + res@DY = nc_file@DY + + if (res@MAP_PROJ .eq. 6) then + res@POLE_LAT = nc_file@POLE_LAT + res@POLE_LON = nc_file@POLE_LON + res@LATINC = (res@DY*360.)/2.0/3.141592653589793/6370000. + res@LONINC = (res@DX*360.)/2.0/3.141592653589793/6370000. + else + res@POLE_LAT = 90.0 + res@POLE_LON = 0.0 + res@LATINC = 0.0 + res@LONINC = 0.0 + end if + + + if(isfilevar(nc_file,"XLAT")) then + if(ISFILE) then + XLAT = nc_file->XLAT(useT,:,:) + XLONG = nc_file->XLONG(useT,:,:) + else + XLAT = file_handle[:]->XLAT + XLONG = file_handle[:]->XLONG + end if + else + if(ISFILE) then + XLAT = nc_file->XLAT_M(useT,:,:) + XLONG = nc_file->XLONG_M(useT,:,:) + else + XLAT = file_handle[:]->XLAT_M + XLONG = file_handle[:]->XLONG_M + end if + end if + + if(dimsizes(dimsizes(XLAT)).eq.2) then +; Rank 2 + res@REF_LAT = XLAT(0,0) + res@REF_LON = XLONG(0,0) + else +; Rank 3 + res@REF_LAT = XLAT(useT,0,0) + res@REF_LON = XLONG(useT,0,0) + end if + res@KNOWNI = 1.0 + res@KNOWNJ = 1.0 + + ; Convert to 1-based indexes for Fortran + new_x = x + 1 + new_y = y + 1 + + loc = wrf_ij_to_ll (new_x,new_y,res) + + if (dimsizes(dimsizes(loc)) .eq. 1) then + loc!0 = "lon_lat" + elseif (dimsizes(dimsizes(loc)) .eq. 2) then + loc!0 = "lon_lat" + loc!1 = "idx" + else ; Not currently supported + loc!0 = "lon_lat" + loc!1 = "domain_idx" + loc!2 = "idx" + end if + + return(loc) + + +end + +;-------------------------------------------------------------------------------- + +undef("wrf_user_vertcross") +function wrf_user_vertcross(var3d:numeric, z_in:numeric, \ + loc_param:numeric, opts:logical ) + +; var3d - 3d field to interpolate (all input fields must be unstaggered) +; z_in - interpolate to this field (either p/z) +; loc_param - an array of 4 values representing the start point and end point +; for the cross section (start_x, start_y, end_x, end_y) OR a single +; point when opt@use_pivot is True representing the pivot point. +; The values can be in grid coordinates or lat/lon coordinates +; (start_x = start_lon, start_y = start_lat, ...). If using +; lat/lon coordinates, then opt@latlon must be True. +; opts - optional arguments +; use_pivot - set to True to indicate that loc_param and angle are used, +; otherwise loc_param is set to 4 values to indicate a start and +; end point +; angle - an angle for vertical plots - 90 represent a WE cross section, +; ignored if use_pivot is False. +; levels - the vertical levels to use in the same units as z_in. Set to +; False to automatically generate the number of levels specified +; by autolevels. +; latlon - set to True if the values in loc_param are latitude and longitude +; values rather than grid values +; file_handle - must be set to a file handle when latlon is True or +; linecoords is True, otherwise this is ignored. +; timeidx - the time index to use for moving nests when latlon is True. Set +; to 0 if the nest is not moving. +; linecoords - set to True to include the latitude and longitude coordinates +; for the cross section line in the output attributes. +; autolevels - set to the desired number of levels when levels are +; selected automatically (default 100). + +begin + + use_pivot = get_res_value(opts, "use_pivot", False) + angle = get_res_value(opts, "angle", 0.0) + levels = get_res_value(opts, "levels", new(1,integer)) + latlon = get_res_value(opts, "latlon", False) + file_handle = get_res_value(opts, "file_handle", 0) + timeidx = get_res_value(opts, "timeidx", 0) + linecoords = get_res_value(opts, "linecoords", False) + nlevels = get_res_value(opts, "autolevels", 100) + + dims = dimsizes(var3d) + nd = dimsizes(dims) + + dimX = dims(nd-1) + dimY = dims(nd-2) + dimZ = dims(nd-3) + + if ( nd .eq. 4 ) then + z = z_in(0,:,:,:) + else + z = z_in + end if + +; Convert latlon to xy coordinates + + if (use_pivot) then + if (latlon) then + opt = True + opt@returnInt = True + opt@useTime = timeidx + ij := wrf_user_ll_to_xy(file_handle, loc_param(0), loc_param(1), opt) + start_x = ij(0) + start_y = ij(1) + else + start_x = loc_param(0) + start_y = loc_param(1) + end if + else + if (latlon) then + opt = True + opt@returnInt = True + opt@useTime = timeidx + ij := wrf_user_ll_to_xy(file_handle, (/ loc_param(0), loc_param(2) /), (/ loc_param(1), loc_param(3) /), opt) + start_x = ij(0,0) + start_y = ij(1,0) + end_x = ij(0,1) + end_y = ij(1,1) + else + start_x = loc_param(0) + start_y = loc_param(1) + end_x = loc_param(2) + end_y = loc_param(3) + end if + end if + +; get the lat/lons along the cross section line if requested + +; set the cross section line coordinates if requested + if (linecoords) then + + latname = "XLAT" + lonname = "XLONG" + if(.not. isfilevar(file_handle,"XLAT")) then + if(isfilevar(file_handle,"XLAT_M")) then + latname = "XLAT_M" + lonname = "XLONG_M" + end if + end if + + latvar = _get_wrf_var(file_handle, latname, timeidx) + lonvar = _get_wrf_var(file_handle, lonname, timeidx) + + if (use_pivot) then + loc := (/start_x, start_y/) + linelats = wrf_user_intrp2d(latvar, loc, angle, False) + linelons = wrf_user_intrp2d(lonvar, loc, angle, False) + else + loc := (/start_x, start_y, end_x, end_y /) + linelats = wrf_user_intrp2d(latvar, loc, angle, True) + linelons = wrf_user_intrp2d(lonvar, loc, angle, True) + end if + + end if + +; set vertical cross section +; Note for wrf_user_set_xy, opt is False when pivot and angle used. + if (use_pivot) then + xy = wrf_user_set_xy( z, start_x, start_y, \ ; assumes 0-based indexing in v6.5.0 + 0.0, 0.0, angle, False ) + + else + xy = wrf_user_set_xy( z, start_x, start_y, \ ; assumes 0-based indexing in v6.5.0 + end_x, end_y, \ + angle, True) + + end if + xp = dimsizes(xy) + + +; first we interp z + var2dz = wrf_interp_2d_xy( z, xy) + +; interp to constant z grid + if (all(ismissing(levels))) then + if(var2dz(0,0) .gt. var2dz(1,0) ) then ; monotonically decreasing coordinate + z_max = floor(max(z)/10)*10 ; bottom value + z_min = ceil(min(z)/10)*10 ; top value + dz = (1.0/nlevels) * (z_max - z_min) + ;nlevels = tointeger( (z_max-z_min)/dz) + z_var2d = new( (/nlevels/), typeof(z)) + z_var2d(0) = z_max + dz = -dz + else + z_max = max(z) + z_min = 0. + dz = (1.0/nlevels) * z_max + ;nlevels = tointeger( z_max/dz ) + z_var2d = new( (/nlevels/), typeof(z)) + z_var2d(0) = z_min + end if + + do i=1, nlevels-1 + z_var2d(i) = z_var2d(0)+i*dz + end do + else + z_var2d = levels + nlevels = dimsizes(z_var2d) + end if + +; interp the variable + if ( dimsizes(dims) .eq. 4 ) then + var2d = new( (/dims(0), nlevels, xp(0)/), typeof(var2dz)) + do it = 0,dims(0)-1 + var2dtmp = wrf_interp_2d_xy( var3d(it,:,:,:), xy) + do i=0,xp(0)-1 + var2d(it,:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d) + end do + end do + var2d!0 = var3d!0 + var2d!1 = "vertical" + var2d!2 = "cross_line_idx" + else + var2d = new( (/nlevels, xp(0)/), typeof(var2dz)) + var2dtmp = wrf_interp_2d_xy( var3d, xy) + do i=0,xp(0)-1 + var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d) + end do + var2d!0 = "vertical" + var2d!1 = "cross_line_idx" + end if + + st_x = tointeger(xy(0,0)) ; + 1 (removed 1-based indexing in 6.5.0) + st_y = tointeger(xy(0,1)) ; + 1 + ed_x = tointeger(xy(xp(0)-1,0)) ; + 1 + ed_y = tointeger(xy(xp(0)-1,1)) ; + 1 + if (.not. use_pivot) then + var2d@Orientation = "Cross-Section: (" + \ + st_x + "," + st_y + ") to (" + \ + ed_x + "," + ed_y + ")" + else + var2d@Orientation = "Cross-Section: (" + \ + st_x + "," + st_y + ") to (" + \ + ed_x + "," + ed_y + ") ; center=(" + \ + start_x + "," + start_y + \ + ") ; angle=" + angle + end if + + if (linecoords) then + var2d@lats = linelats + var2d@lons = linelons + end if + + var2d&vertical = z_var2d + + return(var2d) + +end diff --git a/test/test_filevars.py b/test/test_filevars.py index 1cf3c82..13a04f2 100644 --- a/test/test_filevars.py +++ b/test/test_filevars.py @@ -1,77 +1,79 @@ import unittest as ut -import numpy.testing as nt +import numpy.testing as nt import numpy as np import numpy.ma as ma -import os, sys +import os +import sys import subprocess from wrf import getvar, ALL_TIMES TEST_DIR = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest" TEST_FILENAMES = ["wrfout_d02_2005-08-28_00:00:00", - "wrfout_d02_2005-08-28_12:00:00", - "wrfout_d02_2005-08-29_00:00:00"] + "wrfout_d02_2005-08-28_12:00:00", + "wrfout_d02_2005-08-29_00:00:00"] TEST_FILES = [os.path.join(TEST_DIR, x) for x in TEST_FILENAMES] # Python 3 if sys.version_info > (3,): xrange = range - + class WRFFileVarsTest(ut.TestCase): longMessage = True + def make_test(ncfiles, varname): def test(self): - #import time - #very_start = time.time() - #start = time.time() + # import time + # very_start = time.time() + # start = time.time() t1 = getvar(ncfiles, varname, 0) - - #end = time.time() - #print ("t1: ", start-end) - #start = time.time() + + # end = time.time() + # print("t1: ", start-end) + # start = time.time() t2 = getvar(ncfiles, varname, 0, meta=False) - #end = time.time() - #print ("t2: ", start-end) - #start = time.time() + # end = time.time() + # print("t2: ", start-end) + # start = time.time() t3 = getvar(ncfiles, varname, ALL_TIMES) - #end = time.time() - #print ("t3: ", start-end) - #start = time.time() + # end = time.time() + # print("t3: ", start-end) + # start = time.time() t4 = getvar(ncfiles, varname, ALL_TIMES, meta=False) - #end = time.time() - #print ("t4: ", start-end) - #start = time.time() + # end = time.time() + # print("t4: ", start-end) + # start = time.time() t5 = getvar(ncfiles, varname, ALL_TIMES, method="join") - #end = time.time() - #print ("t5: ", start-end) - #start = time.time() + # end = time.time() + # print("t5: ", start-end) + # start = time.time() t6 = getvar(ncfiles, varname, ALL_TIMES, method="join", meta=False) - #end = time.time() - #print ("t6: ", start-end) - #start = time.time() - - #print ("Total Time: ", (end-start)) + # end = time.time() + # print("t6: ", start-end) + # start = time.time() + + # print("Total Time: ", (end-start)) return test - + if __name__ == "__main__": from netCDF4 import Dataset ncfiles = [Dataset(x) for x in TEST_FILES] - - #import scipy.io - #ncfiles = [scipy.io.netcdf.netcdf_file(x) for x in TEST_FILES] - + + # import scipy.io + # ncfiles = [scipy.io.netcdf.netcdf_file(x) for x in TEST_FILES] + file_vars = ncfiles[0].variables.keys() - + ignore_vars = [] - + for var in file_vars: if var in ignore_vars: continue - + test_func1 = make_test(ncfiles, var) setattr(WRFFileVarsTest, 'test_{0}'.format(var), test_func1) - - ut.main() \ No newline at end of file + + ut.main() diff --git a/test/test_inputs.py b/test/test_inputs.py index 7271681..6f994e4 100644 --- a/test/test_inputs.py +++ b/test/test_inputs.py @@ -1,5 +1,5 @@ import unittest as ut -import numpy.testing as nt +import numpy.testing as nt import numpy as np import numpy.ma as ma import os @@ -17,20 +17,21 @@ TEST_FILES = [os.path.join(IN_DIR, "wrfout_d02_2005-08-28_00:00:00"), os.path.join(IN_DIR, "wrfout_d02_2005-08-28_12:00:00"), os.path.join(IN_DIR, "wrfout_d02_2005-08-29_00:00:00")] + def wrfin_gen(wrf_in): for x in wrf_in: yield x - - + + class wrf_in_iter_class(object): def __init__(self, wrf_in): self._wrf_in = wrf_in self._total = len(wrf_in) self._i = 0 - + def __iter__(self): return self - + def next(self): if self._i >= self._total: raise StopIteration @@ -38,7 +39,7 @@ class wrf_in_iter_class(object): result = self._wrf_in[self._i] self._i += 1 return result - + # Python 3 def __next__(self): return self.next() @@ -48,34 +49,34 @@ def make_test(varname, wrf_in): def test(self): x = getvar(wrf_in, varname, 0) xa = getvar(wrf_in, varname, None) - + return test + def make_interp_test(varname, wrf_in): def test(self): - + # Only testing vinterp since other interpolation just use variables if (varname == "vinterp"): for timeidx in (0, None): eth = getvar(wrf_in, "eth", timeidx=timeidx) - interp_levels = [850,500,5] - field = vinterp(wrf_in, - field=eth, - vert_coord="pressure", - interp_levels=interp_levels, - extrapolate=True, - field_type="theta-e", - timeidx=timeidx, + interp_levels = [850, 500, 5] + field = vinterp(wrf_in, + field=eth, + vert_coord="pressure", + interp_levels=interp_levels, + extrapolate=True, + field_type="theta-e", + timeidx=timeidx, log_p=True) else: pass - - + return test def make_latlon_test(testid, wrf_in): - def test(self): + def test(self): if testid == "xy_out": # Lats/Lons taken from NCL script, just hard-coding for now lats = [-55, -60, -65] @@ -87,42 +88,44 @@ def make_latlon_test(testid, wrf_in): j_s = np.asarray([10, 100, 150], int) ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=0) ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=None) - + return test class WRFVarsTest(ut.TestCase): longMessage = True - + + class WRFInterpTest(ut.TestCase): longMessage = True - + + class WRFLatLonTest(ut.TestCase): longMessage = True + if __name__ == "__main__": - from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, + from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC) omp_set_num_threads(omp_get_num_procs()-1) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) - + ignore_vars = [] # Not testable yet - wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", - "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", - "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", - "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", + wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", + "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", + "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", + "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"] interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"] latlon_tests = ["xy_out", "ll_out"] - for nc_lib in ("netcdf4", "pynio", "scipy"): if nc_lib == "netcdf4": try: from netCDF4 import Dataset except ImportError: - print ("netcdf4-python not installed") + print("netcdf4-python not installed") continue else: test_in = [Dataset(x) for x in TEST_FILES] @@ -130,52 +133,46 @@ if __name__ == "__main__": try: from Nio import open_file except ImportError: - print ("PyNIO not installed") + print("PyNIO not installed") continue else: - test_in = [open_file(x +".nc", "r") for x in TEST_FILES] + test_in = [open_file(x + ".nc", "r") for x in TEST_FILES] elif nc_lib == "scipy": try: from scipy.io.netcdf import netcdf_file except ImportError: - print ("scipy.io.netcdf not installed") + print("scipy.io.netcdf not installed") else: test_in = [netcdf_file(x, mmap=False) for x in TEST_FILES] - + input0 = test_in[0] input1 = test_in input2 = tuple(input1) input3 = wrfin_gen(test_in) input4 = wrf_in_iter_class(test_in) - input5 = {"A" : input1, - "B" : input2} - input6 = {"A" : {"AA" : input1}, - "B" : {"BB" : input2}} - - for i,input in enumerate((input0, input1, input2, + input5 = {"A": input1, + "B": input2} + input6 = {"A": {"AA": input1}, + "B": {"BB": input2}} + + for i, input in enumerate((input0, input1, input2, input3, input4, input5, input6)): for var in wrf_vars: if var in ignore_vars: continue test_func1 = make_test(var, input) - - setattr(WRFVarsTest, "test_{0}_input{1}_{2}".format(nc_lib, - i, var), - test_func1) - + + setattr(WRFVarsTest, "test_{0}_input{1}_{2}".format( + nc_lib, i, var), test_func1) for method in interp_methods: test_interp_func1 = make_interp_test(method, input) - setattr(WRFInterpTest, - "test_{0}_input{1}_{2}".format(nc_lib, - i, method), - test_interp_func1) - + setattr(WRFInterpTest, "test_{0}_input{1}_{2}".format( + nc_lib, i, method), test_interp_func1) + for testid in latlon_tests: test_ll_func = make_latlon_test(testid, input) test_name = "test_{0}_input{1}_{2}".format(nc_lib, i, testid) setattr(WRFLatLonTest, test_name, test_ll_func) - + ut.main() - - \ No newline at end of file diff --git a/test/test_multi_cache.py b/test/test_multi_cache.py index 963fca0..480a70a 100644 --- a/test/test_multi_cache.py +++ b/test/test_multi_cache.py @@ -3,56 +3,60 @@ from wrf.cache import _get_cache from wrf import getvar from netCDF4 import Dataset as nc -#a = nc("/Users/ladwig/Documents/wrf_files/wrf_vortex_single/wrfout_d02_2005-08-28_00:00:00") -#b = nc("/Users/ladwig/Documents/wrf_files/wrf_vortex_single/wrfout_d02_2005-08-28_03:00:00") -a = nc("/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-28_00:00:00") -b = nc("/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-28_12:00:00") -q = {"outoutoutout" : {"outoutout" : {"outout" : {"out1" : {"blah" : [a,b], "blah2" : [a,b]}, "out2" : {"blah" : [a,b], "blah2" : [a,b]} } } } } +a = nc("/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/" + "wrfout_d02_2005-08-28_00:00:00") +b = nc("/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/" + "wrfout_d02_2005-08-28_12:00:00") +q = {"outoutoutout": + {"outoutout": + {"outout": + {"out1": {"blah": [a, b], "blah2": [a, b]}, + "out2": {"blah": [a, b], "blah2": [a, b]}}}}} t1 = time.time() c = getvar(q, "rh", method="cat", timeidx=None, squeeze=True) t2 = time.time() -print (c) -print ("time taken: {}".format((t2-t1)*1000.)) +print(c) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="cat", timeidx=None, squeeze=False) t2 = time.time() -print (c) -print ("time taken: {}".format((t2-t1)*1000.)) +print(c) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="cat", timeidx=1, squeeze=True) t2 = time.time() -print (c) -print ("time taken: {}".format((t2-t1)*1000.)) +print(c) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="cat", timeidx=1, squeeze=False) t2 = time.time() print(c) -print ("time taken: {}".format((t2-t1)*1000.)) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="join", timeidx=None, squeeze=True) t2 = time.time() -print (c) -print ("time taken: {}".format((t2-t1)*1000.)) +print(c) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="join", timeidx=None, squeeze=False) t2 = time.time() print(c) -print ("time taken: {}".format((t2-t1)*1000.)) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="join", timeidx=1, squeeze=True) t2 = time.time() -print (c) -print ("time taken: {}".format((t2-t1)*1000.)) +print(c) +print("time taken: {}".format((t2-t1)*1000.)) t1 = time.time() c = getvar(q, "rh", method="join", timeidx=1, squeeze=False) t2 = time.time() -print (c) -print ("time taken: {}".format((t2-t1)*1000.)) +print(c) +print("time taken: {}".format((t2-t1)*1000.)) diff --git a/test/test_omp.py b/test/test_omp.py index 5c621a9..718ab2f 100644 --- a/test/test_omp.py +++ b/test/test_omp.py @@ -1,16 +1,16 @@ -from __future__ import (absolute_import, division, print_function, +from __future__ import (absolute_import, division, print_function, unicode_literals) import unittest as ut -import numpy.testing as nt +import numpy.testing as nt -from wrf import (omp_set_num_threads, omp_get_num_threads, +from wrf import (omp_set_num_threads, omp_get_num_threads, omp_get_max_threads, omp_get_thread_num, - omp_get_num_procs, omp_in_parallel, + omp_get_num_procs, omp_in_parallel, omp_set_dynamic, omp_get_dynamic, omp_set_nested, - omp_get_nested, omp_set_schedule, - omp_get_schedule, omp_get_thread_limit, - omp_set_max_active_levels, + omp_get_nested, omp_set_schedule, + omp_get_schedule, omp_get_thread_limit, + omp_set_max_active_levels, omp_get_max_active_levels, omp_get_level, omp_get_ancestor_thread_num, omp_get_team_size, omp_get_active_level, omp_in_final, @@ -20,103 +20,102 @@ from wrf import (omp_set_num_threads, omp_get_num_threads, omp_unset_lock, omp_unset_nest_lock, omp_test_lock, omp_test_nest_lock, omp_get_wtime, omp_get_wtick, - OMP_SCHED_STATIC, OMP_SCHED_DYNAMIC, - OMP_SCHED_GUIDED, OMP_SCHED_AUTO) - + OMP_SCHED_STATIC, OMP_SCHED_DYNAMIC, + OMP_SCHED_GUIDED, OMP_SCHED_AUTO) + class OmpTest(ut.TestCase): longMessage = True - + def test_locks(self): - l = omp_init_lock() - omp_set_lock(l) - omp_unset_lock(l) - omp_test_lock(l) - omp_destroy_lock(l) - + lk = omp_init_lock() + omp_set_lock(lk) + omp_unset_lock(lk) + omp_test_lock(lk) + omp_destroy_lock(lk) + nl = omp_init_nest_lock() omp_set_nest_lock(nl) omp_unset_nest_lock(nl) omp_test_nest_lock(nl) omp_destroy_nest_lock(nl) - - + def test_thread_set(self): omp_set_num_threads(4) max_threads = omp_get_max_threads() self.assertEqual(max_threads, 4) - + num_threads = omp_get_num_threads() - self.assertEqual(num_threads, 1) # Always 1 outside of parallel region - + # Always 1 outside of parallel region + self.assertEqual(num_threads, 1) + thread_num = omp_get_thread_num() - self.assertEqual(thread_num, 0) # Always 0 outside of parallel region + # Always 0 outside of parallel region + self.assertEqual(thread_num, 0) num_procs = omp_get_num_procs() in_parallel = omp_in_parallel() - self.assertFalse(in_parallel) # Always False outside of parallel region - + # Always False outside of parallel region + self.assertFalse(in_parallel) + limit = omp_get_thread_limit() - - + def test_dynamic(self): omp_set_dynamic(True) dynamic = omp_get_dynamic() self.assertTrue(dynamic) - + omp_set_dynamic(False) dynamic = omp_get_dynamic() self.assertFalse(dynamic) - + def test_nested(self): omp_set_nested(True) nested = omp_get_nested() self.assertTrue(nested) - + omp_set_nested(False) nested = omp_get_nested() self.assertFalse(nested) - - + def test_schedule(self): omp_set_schedule(OMP_SCHED_STATIC, 100000) kind, modifier = omp_get_schedule() self.assertEqual(kind, OMP_SCHED_STATIC) self.assertEqual(modifier, 100000) - + omp_set_schedule(OMP_SCHED_DYNAMIC, 10000) kind, modifier = omp_get_schedule() self.assertEqual(kind, OMP_SCHED_DYNAMIC) self.assertEqual(modifier, 10000) - - omp_set_schedule(OMP_SCHED_GUIDED, 100) + + omp_set_schedule(OMP_SCHED_GUIDED, 100) kind, modifier = omp_get_schedule() self.assertEqual(kind, OMP_SCHED_GUIDED) self.assertEqual(modifier, 100) - - omp_set_schedule(OMP_SCHED_AUTO, 10) + + omp_set_schedule(OMP_SCHED_AUTO, 10) kind, modifier = omp_get_schedule() self.assertEqual(kind, OMP_SCHED_AUTO) - self.assertNotEqual(modifier, 10) # The modifier argument is ignored, - # so it will be set to the previous - # value of 100. - - + # The modifier argument is ignored, + # so it will be set to the previous + # value of 100. + self.assertNotEqual(modifier, 10) + def test_team_level(self): omp_set_max_active_levels(10) active_levels = omp_get_max_active_levels() self.assertEqual(active_levels, 10) - + level = omp_get_level() ancestor_thread = omp_get_ancestor_thread_num(level) team_size = omp_get_team_size(level) active_level = omp_get_active_level() in_final = omp_in_final() - - + def test_time(self): wtime = omp_get_wtime() wtick = omp_get_wtick() + if __name__ == "__main__": ut.main() - \ No newline at end of file diff --git a/test/test_proj_params.py b/test/test_proj_params.py index 77a0ec7..460fae5 100644 --- a/test/test_proj_params.py +++ b/test/test_proj_params.py @@ -1,314 +1,302 @@ import unittest as ut -import numpy.testing as nt +import numpy.testing as nt import numpy as np import numpy.ma as ma -import os, sys +import os +import 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,): + 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, + 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 - + + 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, + 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, + 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.) - + 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.) - + 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, + + 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, + 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.) - + 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, + + 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.) - + 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, + + 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, + 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, + 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, + 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.) - + 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.*" ): + 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.*" ): + + 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, + + 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, + + 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 - + + # 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, + 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, + 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.) - + 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.) - + 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, + + 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, + 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.) - + 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, + + 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.) - + 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, + + 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, + 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, + 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, + 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.) - + 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/test_units.py b/test/test_units.py index 40c07eb..f1cab6f 100644 --- a/test/test_units.py +++ b/test/test_units.py @@ -1,7 +1,7 @@ import sys import unittest as ut -import numpy.testing as nt +import numpy.testing as nt import numpy as np import numpy.ma as ma from netCDF4 import Dataset as nc @@ -13,44 +13,43 @@ TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00" # Python 3 if sys.version_info > (3,): xrange = range - + class TestUnits(ut.TestCase): longMessage = True - + def test_temp_units(self): wrfnc = nc(TEST_FILE) - + for units in ("K", "degC", "degF"): var = getvar(wrfnc, "temp", units=units) - + self.assertEqual(var.attrs["units"], units) - + def test_wind_units(self): wrfnc = nc(TEST_FILE) - + for units in ("m s-1", "kt", "mi h-1", "km h-1", "ft s-1"): var = getvar(wrfnc, "uvmet", units=units) - + self.assertEqual(var.attrs["units"], units) - + def test_pres_units(self): wrfnc = nc(TEST_FILE) - + for units in ("Pa", "hPa", "mb", "torr", "mmHg", "atm"): var = getvar(wrfnc, "slp", units=units) - + self.assertEqual(var.attrs["units"], units) - + def test_height_units(self): wrfnc = nc(TEST_FILE) - + for units in ("m", "km", "dam", "ft", "mi"): var = getvar(wrfnc, "z", units=units) - + self.assertEqual(var.attrs["units"], units) - - + if __name__ == "__main__": - ut.main() \ No newline at end of file + ut.main() diff --git a/test/utests.py b/test/utests.py index 08a1087..eb95167 100644 --- a/test/utests.py +++ b/test/utests.py @@ -35,7 +35,7 @@ def setUpModule(): os.environ["OMP_NUM_THREADS"] = "4" this_path = os.path.realpath(__file__) - ncl_script = os.path.join(os.path.dirname(this_path), + ncl_script = os.path.join(os.path.dirname(this_path), "ncl", "ncl_get_var.ncl") for dir, outfile in zip(DIRS, REF_NC_FILES):