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):