From 32dffa285fa26393ec21ad0c7ddc4c1e9dbfa8a1 Mon Sep 17 00:00:00 2001 From: Nick Guy Date: Wed, 21 Feb 2018 00:04:13 -0800 Subject: [PATCH 1/3] PEP8 and initial travis CI --- .travis.yml | 40 ++++ README.md | 4 +- pyradarmet/BeamBlock.py | 431 +++++++++++++++++++--------------------- pyradarmet/__init__.py | 2 +- pyradarmet/doppler.py | 14 +- pyradarmet/geometry.py | 12 +- pyradarmet/system.py | 37 ++-- pyradarmet/variables.py | 12 +- pyradarmet/zdrcal.py | 297 +++++++++++++-------------- scripts/cal_zdr | 142 +++++++------ scripts/listfile_raw.py | 65 +++--- scripts/zdrcal_3panel | 96 ++++----- setup.py | 80 ++++---- 13 files changed, 632 insertions(+), 600 deletions(-) create mode 100644 .travis.yml diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..2fd4d04 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,40 @@ +sudo: false # Use container-based infrastructure +language: python + +matrix: + include: + - python: 2.7 + env: + - PYTHON_VERSION="2.7" + - NOSE_ARGS="-v --with-cov --cov pyradarmet pyradarmet" + - COVERALLS="true" + - python: 2.7 + env: + - PYTHON_VERSION="2.7" + - NOSE_ARGS="-v --with-cov --cov pyradarmet --exe pyradarmet" + - FROM_RECIPE="true" + - DOC_BUILD="true" + - python: 3.4 + env: + - PYTHON_VERSION="3.4" + - NOSE_ARGS="-v --with-cov --cov pyradarmet pyradarmet" + - python: 3.5 + env: + - PYTHON_VERSION="3.5" + - NOSE_ARGS="-v --with-cov --cov pyradarmet pyradarmet" + - python: 3.6 + env: + - PYTHON_VERSION="3.6" + - NOSE_ARGS="-v --with-cov --cov pyradarmet pyradarmet" +before_install: +- pip install pytest pytest-cov +- pip install coveralls +install: source continuous_integration/install.sh +script: eval xvfb-run nosetests $NOSE_ARGS +after_success: + # Ignore coveralls failures as the coveralls server is not very reliable + # but we don't want travis to report a failure in the github UI just + # because the coverage report failed to be published. + - if [[ "$COVERALLS" == "true" ]]; then coveralls || echo "failed"; fi + # Build docs if requested +- if [[ "$DOC_BUILD" == "true" ]]; then cd $TRAVIS_BUILD_DIR; source continuous_integration/build_docs.sh; fi \ No newline at end of file diff --git a/README.md b/README.md index 991edb8..ffc1708 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,11 @@ PyRadarMet -=============== +========== Python Fundamental Calculations in Radar Meteorology package notes Originally Created: 5 February 2014 ## Author -Nick Guy - nick.guy@uwyo.edu +Nick Guy - nguysur@gmail.com Special thanks to Timothy Lang and Kai Muehlbauer for the insights and contributions. diff --git a/pyradarmet/BeamBlock.py b/pyradarmet/BeamBlock.py index a4c7339..6e686f9 100644 --- a/pyradarmet/BeamBlock.py +++ b/pyradarmet/BeamBlock.py @@ -5,8 +5,8 @@ Author:: Nick Guy - OU CIMMS/Univ of Miami -This program was ported from code written in IDL used at Colorado State University. -It is believed that the original code at CSU was written by Steven Nesbitt. +This program was ported from IDL code used at Colorado State University. +It is believed to be originally written at CSU by Steven Nesbitt. Timothy Lang also contributed to the development of this program. This program is not particularly fast as it reads in large DEM files. @@ -49,16 +49,14 @@ PAN3_AX3 = [.10, .08, .80, .35] # Set up the data directory for DEM files -#DEM_DATADIR = os.sep.join([os.path.dirname(__file__), 'data']) -#print DEM_DATADIR +# DEM_DATADIR = os.sep.join([os.path.dirname(__file__), 'data']) +# print DEM_DATADIR + class BeamBlock(object): """ A class instance to calculalate beam blockage of weather radar """ - -################################################### - def __init__(self, dem_filename=None, verbose=False, dem_hdr_file=None, upper_left_x=None, upper_left_y=None, radar_lon=None, radar_lat=None, radar_alt=None, @@ -70,19 +68,19 @@ def __init__(self, dem_filename=None, verbose=False, """ If initialized with a filename (incl. path), will call read_bin_netcdf() to populate the class instance. - + If not, it simply instances the class but does not populate its attributes. - + Parameters:: ---------- dem_filename : str Full path and filename of Digital Elevation Model file to use. - If not supplied, a file will be selected, along with the + If not supplied, a file will be selected, along with the appropriate dem_hdr_file, based on radar_lon and radar_lat dem_hdr_file : str - Full path and filename of Digital Elevation Model header file to use. - This file will find upper_left_x and upper_left_y. + Full path and filename of Digital Elevation Model header file. + This file will find upper_left_x and upper_left_y. If not set, upper_left_x and upper_left_y must be set by user. upper_left_x : float Coordinate of upper left X position of file [decimal degrees] @@ -91,7 +89,7 @@ def __init__(self, dem_filename=None, verbose=False, Coordinate of upper left Y position of file [decimal degrees] verbose: boolean Set to True for text output. Useful for debugging. - + radar_lon : float Radar longitude [decimal degrees] radar_lat : float @@ -127,29 +125,29 @@ def __init__(self, dem_filename=None, verbose=False, self.lon_box_size = lon_box_size self.lat_box_size = lat_box_size self.ralt = radar_alt - + # Set the directory where DEM data can be found try: os.environ["GTOPO_DATA"] self.demdir = os.path.join(os.environ["GTOPO_DATA"], 'data') except KeyError: - print "No GTOPO_DATA environmental variable found, "\ - "assuming data embedded in PyRadarMet package" + print("No GTOPO_DATA environmental variable found, " + "assuming data embedded in PyRadarMet package") self.demdir = os.sep.join([os.path.dirname(__file__), 'data']) - + # Set radar latitude and longitude if (radar_lon is None) or (radar_lat is None): print "Radar longitude, latitude, and altitude need to be defined!" exit() else: - self.rlon, self.rlat= radar_lon, radar_lat - + self.rlon, self.rlat = radar_lon, radar_lat + # Set the height of the radar antenna above ground if radar_antenna_height is None: self.rheight = 0. else: self.rheight = radar_antenna_height - + # Set elevation angle (degrees) [e.g. 0.8, 1.3, 1.8, 2.5] if elev_angle is None: self.E = 0.8 @@ -185,24 +183,22 @@ def __init__(self, dem_filename=None, verbose=False, self.nrng = 601 else: self.nrng = num_rng_pts - + # Set the number of azimuthal points if num_az_pts is None: self.naz = 361 else: self.naz = num_az_pts - + # Check the DEM file self._check_dem_file() # self._grab_dem_file() - + # Read the DEM file self.read_dem() - + # Process inputs to calculate beam properties self.calc_beam_blockage() - - ############## def help(self): @@ -220,10 +216,6 @@ def help(self): print 'diag(), get_comp(), plot_vert(), plot_horiz(), three_panel_plot()' print 'subsection(), write_mosaic_binary(), output_composite()' _method_footer_printout() - -########################## -# Read data file methods # -########################## def read_dem(self): '''Read a DEM file if passed or search for correct file @@ -240,8 +232,9 @@ def read_gt30_dem(self): The following code was adapted from an example at: http://stevendkay.wordpress.com/tag/dem/ ''' - if isinstance(self.dem_hdr_file, str) != False: - # wrap header in dict, because the different elements aren't necessarily in the same row + if isinstance(self.dem_hdr_file, str) is not False: + # wrap header in dict, because the different elements aren't + # necessarily in the same row hdr = dict(np.loadtxt(self.dem_hdr_file, dtype='S20')) byteorder = hdr['BYTEORDER'] pixeltype = hdr['PIXELTYPE'] @@ -263,7 +256,7 @@ def read_gt30_dem(self): missing = -9999. xdim = 0.00833333333333 ydim = 0.00833333333333 - + # Extract GTOPO30 BIL file fi = open(self.dem_file, "rb") databin = fi.read() @@ -280,72 +273,67 @@ def read_gt30_dem(self): s = ">%d%s" % (nrows * ncols, format,) else: s = "<%d%s" % (nrows * ncols, format,) - z = struct.unpack(s, databin) - + topo = np.zeros((nrows, ncols)) for r in range(0, nrows): for c in range(0, ncols): elevation = z[((ncols) * r) + c] if (elevation == 65535 or elevation < 0 or elevation > 20000): - # may not be needed depending on format, and the "magic number" - # value used for 'void' or missing data - elevation = 0.#np.nan + # may not be needed depending on format, and the + # "magic number" value used for 'void' or missing data + elevation = 0. topo[r][c] = float(elevation) lat = upper_left_y - ydim * np.arange(nrows) lon = upper_left_x + xdim * np.arange(ncols) - + # Find indices for a subset of DEM data londel = self.lon_box_size / 2. latdel = self.lat_box_size / 2. - lonInd = np.where((lon >= (self.rlon - londel)) & \ + lonInd = np.where((lon >= (self.rlon - londel)) & (lon <= (self.rlon + londel))) - latInd = np.where((lat >= (self.rlat - latdel)) & \ + latInd = np.where((lat >= (self.rlat - latdel)) & (lat <= (self.rlat + latdel))) lonSub = lon[np.min(lonInd):np.max(lonInd)] latSub = lat[np.min(latInd): np.max(latInd)] self.lon, self.lat = np.meshgrid(lonSub, latSub) - self.topo = topo[np.min(latInd):np.max(latInd), np.min(lonInd): np.max(lonInd)] - + self.topo = topo[np.min(latInd):np.max(latInd), + np.min(lonInd):np.max(lonInd)] + self.topo = np.ma.masked_outside(self.topo, 0, 20000) self.topo = np.ma.masked_invalid(self.topo) - + # Find the elevation of the radar if self.ralt is None: radar_lon_index = self._get_array_index(self.rlon, self.lon[0, :]) radar_lat_index = self._get_array_index(self.rlat, self.lat[:, 0]) - - radar_alt_map = scim.map_coordinates(self.topo, - np.vstack((radar_lat_index, radar_lon_index)), - prefilter=False) - self.ralt = self.rheight + radar_alt_map -################################################### -# Calculation methods # -################################################### + radar_alt_map = scim.map_coordinates( + self.topo, np.vstack((radar_lat_index, radar_lon_index)), + prefilter=False) + self.ralt = self.rheight + radar_alt_map def calc_beam_blockage(self): '''Calculate the properties of the beam for calculation''' # Initialize arrays self.inX = np.empty(self.nrng) self.inY = np.empty(self.nrng) - self.rng_lon = np.empty([self.naz, self.nrng]) self.rng_lat = np.empty([self.naz, self.nrng]) self.terr = np.empty([self.naz, self.nrng]) - + # Points along each ray to calculate beam propagation # e.g. From 0 to 150 km, with 0.25 km bin spacing self.rng = np.linspace(0, self.range*1000., self.nrng) - + # Calculate the effective radius self.Reff = r_effective() # Calculate the height of center of beam - self.h = ray_height(self.rng, self.E, - self.ralt, reff=self.Reff) + self.h = ray_height( + self.rng, self.E, self.ralt, reff=self.Reff) # Calculate the beam radius self.a = half_power_radius(self.rng, self.BW) @@ -355,44 +343,38 @@ def calc_beam_blockage(self): # Set up arrays for calculations self.phis = np.linspace(0, 360, self.naz) - + # Initialize the beam blockage arrays with 0 self.PBB = np.zeros((self.naz, self.nrng)) self.CBB = np.zeros((self.naz, self.nrng)) for jj, az in enumerate(self.phis): # Calculate the Lons/Lats along the rng_gnd - self.rng_lat[jj, :] = np.degrees(np.arcsin(np.sin(np.radians(self.rlat)) * \ - np.cos(self.rng_gnd / RE) +\ - np.cos(np.radians(self.rlat)) * \ - np.sin(self.rng_gnd / RE) *\ - np.cos(np.radians(az))\ - )\ - ) + self.rng_lat[jj, :] = np.degrees( + np.arcsin(np.sin(np.radians(self.rlat)) * + np.cos(self.rng_gnd / RE) + + np.cos(np.radians(self.rlat)) * + np.sin(self.rng_gnd / RE) * + np.cos(np.radians(az)))) self.rng_lon[jj, :] = self.rlon + \ - np.degrees(np.arctan2(np.sin(np.radians(az)) * \ - np.sin(self.rng_gnd / RE) *\ - np.cos(np.radians(self.rlat)), \ - np.cos(self.rng_gnd / RE) - \ - np.sin(np.radians(self.rlat)) * \ - np.sin(np.radians(self.rng_lat[jj, :]))\ - )\ - ) - + np.degrees(np.arctan2(np.sin(np.radians(az)) * + np.sin(self.rng_gnd / RE) * + np.cos(np.radians(self.rlat)), + np.cos(self.rng_gnd / RE) - + np.sin(np.radians(self.rlat)) * + np.sin(np.radians(self.rng_lat[jj, :])))) + # Find the indices for interpolation for ii, junk in enumerate(self.inX): self.inX[ii] = self._get_array_index(self.rng_lon[jj, ii], self.lon[0, :]) self.inY[ii] = self._get_array_index(self.rng_lat[jj, ii], self.lat[:, 0]) - # Interpolate terrain heights to the radar grid - self.terr[jj, :] = scim.map_coordinates(self.topo, - np.vstack((self.inY, self.inX)), - prefilter=False) + self.terr[jj, :] = scim.map_coordinates( + self.topo, np.vstack((self.inY, self.inX)), prefilter=False) # Calculate PBB along range - self.PBB[jj, :] = beam_block_frac(self.terr[jj, :], - self.h, self.a) + self.PBB[jj, :] = beam_block_frac(self.terr[jj, :], self.h, self.a) self.PBB = np.ma.masked_invalid(self.PBB) @@ -400,13 +382,10 @@ def calc_beam_blockage(self): if ii == 0: self.CBB[:, ii] = self.PBB[:, ii] else: - self.CBB[:,ii] = np.fmax([self.PBB[:, ii]], [self.CBB[:, ii-1]]) - - self.terr = np.ma.masked_less(self.terr, -1000.) + self.CBB[:, ii] = np.fmax([self.PBB[:, ii]], + [self.CBB[:, ii-1]]) -############### -# Get methods # -############### + self.terr = np.ma.masked_less(self.terr, -1000.) def _get_array_index(self, value, array): '''Calculate the exact index position within latitude array''' @@ -416,24 +395,19 @@ def _get_array_index(self, value, array): pos = np.absolute(value - array[0]) / dp return pos -#################### -# DEM file methods # -#################### - def _check_dem_file(self): '''Check if dem file is passed and set if not''' - ###Add ability to stitch 2 files together self.rlon - londel + # Add ability to stitch 2 files together self.rlon - londel # Set the radar min and max longitude to plot, arbitrary chosen minlat, maxlat = self.rlat - 4, self.rlat + 4 minlon, maxlon = self.rlon - 4, self.rlon + 4 # get radar bounding box corner coords radar_corners = np.array([[minlon, minlat], [minlon, maxlat], - [maxlon, maxlat], [maxlon, minlat]]) + [maxlon, maxlat], [maxlon, minlat]]) # file list which will take the found hdr files flist = [] - # iterate over dem files, besides antarcps for dem_hdr_file in sorted(glob(os.path.join(self.demdir, '[e,w]*.hdr'))): hdr = dict(np.loadtxt(dem_hdr_file, dtype='S20')) @@ -449,8 +423,10 @@ def _check_dem_file(self): lower_right_y = upper_left_y - nrows * ydim # construct dem bounding box - dem_bbox = [[upper_left_x, lower_right_y], [upper_left_x, upper_left_y], - [lower_right_y, upper_left_y], [lower_right_x, lower_right_y]] + dem_bbox = [[upper_left_x, lower_right_y], + [upper_left_x, upper_left_y], + [lower_right_y, upper_left_y], + [lower_right_x, lower_right_y]] # make matplotlib Path from dem bounding box dem_bbox_path = Path(dem_bbox) @@ -466,8 +442,9 @@ def _check_dem_file(self): break # construct gdalbuildvirt command - command = ['gdalbuildvrt', '-overwrite', '-te', str(minlon), str(minlat), - str(maxlon), str(maxlat), os.path.join(self.demdir, 'interim.vrt')] + command = ['gdalbuildvrt', '-overwrite', '-te', str(minlon), + str(minlat), str(maxlon), str(maxlat), + os.path.join(self.demdir, 'interim.vrt')] for f in flist: command.append(f) @@ -486,128 +463,125 @@ def _check_dem_file(self): if self.dem_file is None: if (self.rlon >= -180.) and (self.rlon < 180.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'antarcps' elif (self.rlon >= 20.) and (self.rlon < 60.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'e020n40' elif (self.rlon >= 20.) and (self.rlon < 60.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'e020n90' elif (self.rlon >= 20.) and (self.rlon < 60.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'e020s10' - + elif (self.rlon >= 60.) and (self.rlon < 100.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'e060n40' elif (self.rlon >= 60.) and (self.rlon < 100.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'e060n90' elif (self.rlon >= 60.) and (self.rlon < 100.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'e060s10' elif (self.rlon >= 60.) and (self.rlon < 120.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'e060s60' - + elif (self.rlon >= 100.) and (self.rlon < 140.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'e100n40' elif (self.rlon >= 100.) and (self.rlon < 140.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'e100n90' elif (self.rlon >= 100.) and (self.rlon < 140.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'e100s10' elif (self.rlon >= 120.) and (self.rlon < 180.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'e120s60' - + elif (self.rlon >= 140.) and (self.rlon < 180.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'e140n40' elif (self.rlon >= 140.) and (self.rlon < 180.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'e140n90' elif (self.rlon >= 140.) and (self.rlon < 180.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'e140s10' - + elif (self.rlon >= 0.) and (self.rlon < 60.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'w000s60' - + elif (self.rlon >= -20.) and (self.rlon < 20.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'w020n40' elif (self.rlon >= -20.) and (self.rlon < 20.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'w020n90' elif (self.rlon >= -20.) and (self.rlon < 20.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'w020s10' - + elif (self.rlon >= -60.) and (self.rlon < -20.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'w060n40' elif (self.rlon >= -60.) and (self.rlon < -20.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'w060n90' elif (self.rlon >= -60.) and (self.rlon < -20.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'w060s10' elif (self.rlon >= -60.) and (self.rlon < 0.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'w060s60' - + elif (self.rlon >= -100.) and (self.rlon < -60.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'w100n40' elif (self.rlon >= -100.) and (self.rlon < -60.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'w100n90' elif (self.rlon >= -100.) and (self.rlon < -60.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'w100s10' elif (self.rlon >= -120.) and (self.rlon < -60.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'w120s60' - + elif (self.rlon >= -140.) and (self.rlon < -100.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'w140n40' elif (self.rlon >= -140.) and (self.rlon < -100.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'w140n90' elif (self.rlon >= -140.) and (self.rlon < -100.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'w140s10' - + elif (self.rlon >= -180.) and (self.rlon < -140.) and \ - (self.rlat >= -10.) and (self.rlat < 40.): + (self.rlat >= -10.) and (self.rlat < 40.): dem_filename = 'w180n40' elif (self.rlon >= -180.) and (self.rlon < -140.) and \ - (self.rlat >= 40.) and (self.rlat < 90.): + (self.rlat >= 40.) and (self.rlat < 90.): dem_filename = 'w180n90' elif (self.rlon >= -180.) and (self.rlon < -140.) and \ - (self.rlat >= -60.) and (self.rlat < -10.): + (self.rlat >= -60.) and (self.rlat < -10.): dem_filename = 'w180s10' elif (self.rlon >= -180.) and (self.rlon < -120.) and \ - (self.rlat >= -90.) and (self.rlat < -60.): + (self.rlat >= -90.) and (self.rlat < -60.): dem_filename = 'w180s60' - + self.dem_file = os.path.join(self.demdir, (dem_filename + ".dem")) - self.dem_hdr_file = os.path.join(self.demdir, (dem_filename + ".hdr")) - - # Need an alternative method here to find points close to boundary - dem_filename2 = None + self.dem_hdr_file = os.path.join( + self.demdir, (dem_filename + ".hdr")) + + # Need an alternative method here to find points close to boundary + dem_filename2 = None if dem_filename2 is not None: print "MAY NEED A SECOND DEM FILE" - -################ -# Plot methods # -################ def plot_3panel(self, beam_blockage='complete', lon_plot_range=2., lat_plot_range=2., @@ -620,8 +594,8 @@ def plot_3panel(self, beam_blockage='complete', saveFig=False): """ Create a 3-panel plot characterization of beam blockage - - Parameters:: + + Parameters ---------- beam_blockage : str 'complete' or 'partial' chooses which field to plot @@ -635,9 +609,11 @@ def plot_3panel(self, beam_blockage='complete', beam_blockage_cmap : str Matplotlib colormap to use for beam blockage map elev_min : float - Minumum elevation to display on beam propagation and terrain height [meters] + Minumum elevation to display on beam propagation and + terrain height [meters] elev_max : float - Maximum elevation to display on beam propagation and terrain height [meters] + Maximum elevation to display on beam propagation and + terrain height [meters] range_rings : float A list with location of range rings - e.g. [50., 100.] ht_shade_min : float @@ -651,36 +627,38 @@ def plot_3panel(self, beam_blockage='complete', saveFig : boolean True to save figure as output file, False to show """ - # Set the radar min and max longitude to plot - self.minlat, self.maxlat = self.rlat - lat_plot_range, self.rlat + lat_plot_range - self.minlon, self.maxlon = self.rlon - lon_plot_range, self.rlon + lon_plot_range + self.minlat = self.rlat - lat_plot_range + self.maxlat = self.rlat + lat_plot_range + self.minlon = self.rlon - lon_plot_range + self.maxlon = self.rlon + lon_plot_range if terrain_cmap is None: terrain_cmap = 'BrBG_r' if beam_blockage_cmap is None: beam_blockage_cmap = 'PuRd' self.terr_cmap, self.bb_cmap = terrain_cmap, beam_blockage_cmap - - fig = plt.figure(figsize=(10,8)) + + fig = plt.figure(figsize=(10, 8)) ax1 = fig.add_axes(PAN3_AX1) ax2 = fig.add_axes(PAN3_AX2) ax3 = fig.add_axes(PAN3_AX3) - - self.draw_terrain_height_map(fig, ax1, vmin=ht_shade_min, vmax=ht_shade_max, - lat_spacing=lat_spacing, lon_spacing=lon_spacing) + + self.draw_terrain_height_map( + fig, ax1, vmin=ht_shade_min, vmax=ht_shade_max, + lat_spacing=lat_spacing, lon_spacing=lon_spacing) if beam_blockage == 'partial': self.draw_bb_map(fig, ax2, BB=self.PBB, range_rings=range_rings, - lat_spacing=lat_spacing, lon_spacing=lon_spacing) + lat_spacing=lat_spacing, lon_spacing=lon_spacing) elif beam_blockage == 'complete': self.draw_bb_map(fig, ax2, range_rings=range_rings, - lat_spacing=lat_spacing, lon_spacing=lon_spacing) + lat_spacing=lat_spacing, lon_spacing=lon_spacing) self.draw_beam_terrain_profile(fig, ax3, ymin=elev_min, ymax=elev_max) - + if saveFig: plt.savefig('BBmap.png', format='png') else: plt.show() - + def draw_terrain_height_map(self, fig, ax, vmin=None, vmax=None, lat_spacing=None, lon_spacing=None): '''Draw the terrain heights''' @@ -692,48 +670,51 @@ def draw_terrain_height_map(self, fig, ax, vmin=None, vmax=None, topomax = 3. else: topomax = vmax/1000. - + if lat_spacing is None: lat_spacing = 1. if lon_spacing is None: lon_spacing = 1. - - bm1 = Basemap(projection='cea', resolution='l', area_thresh = 10000., - llcrnrlon=self.minlon, urcrnrlon=self.maxlon, - llcrnrlat=self.minlat, urcrnrlat=self.maxlat, - ax=ax) - ax.set_title('Terrain within %02d km of Radar (km)'%(self.range), fontdict=TITLEDICT) - bm1.drawmeridians(np.arange(self.minlon, self.maxlon, lon_spacing), labels=[1,0,0,1]) - bm1.drawparallels(np.arange(self.minlat, self.maxlat, lat_spacing), labels=[1,0,0,1]) + + bm1 = Basemap(projection='cea', resolution='l', area_thresh=10000., + llcrnrlon=self.minlon, urcrnrlon=self.maxlon, + llcrnrlat=self.minlat, urcrnrlat=self.maxlat, ax=ax) + ax.set_title('Terrain within %02d km of Radar (km)' % (self.range), + fontdict=TITLEDICT) + bm1.drawmeridians(np.arange(self.minlon, self.maxlon, lon_spacing), + labels=[1, 0, 0, 1]) + bm1.drawparallels(np.arange(self.minlat, self.maxlat, lat_spacing), + labels=[1, 0, 0, 1]) bm1.drawcountries() bm1.drawcoastlines() bm1.drawrivers() - + xbm1, ybm1 = bm1(self.lon, self.lat) - Htmap = bm1.pcolormesh(xbm1, ybm1, self.topo/1000., vmin=topomin, vmax=topomax, - cmap=self.terr_cmap) + Htmap = bm1.pcolormesh(xbm1, ybm1, self.topo/1000., vmin=topomin, + vmax=topomax, cmap=self.terr_cmap) bm1.plot(self.rlon, self.rlat, 'rD', latlon=True) - #plot_range_ring(50., bm=bm1, color='w') + # plot_range_ring(50., bm=bm1, color='w') fig.colorbar(Htmap, ax=ax) - + def draw_bb_map(self, fig, ax, BB=None, range_rings=None, - lat_spacing=None, lon_spacing=None): + lat_spacing=None, lon_spacing=None): '''Draw the Beam Blockage''' if BB is None: BB = self.CBB - + if lat_spacing is None: lat_spacing = 1. if lon_spacing is None: lon_spacing = 1. - - bm2 = Basemap(projection='cea', resolution='l', area_thresh = 10000., - llcrnrlon=self.minlon, urcrnrlon=self.maxlon, - llcrnrlat=self.minlat, urcrnrlat=self.maxlat, - ax=ax) + + bm2 = Basemap(projection='cea', resolution='l', area_thresh=10000., + llcrnrlon=self.minlon, urcrnrlon=self.maxlon, + llcrnrlat=self.minlat, urcrnrlat=self.maxlat, ax=ax) ax.set_title('Beam-blockage fraction', fontdict=TITLEDICT) - bm2.drawmeridians(np.arange(self.minlon, self.maxlon, lon_spacing), labels=[1,0,0,1]) - bm2.drawparallels(np.arange(self.minlat, self.maxlat, lat_spacing), labels=[1,0,0,1]) + bm2.drawmeridians(np.arange(self.minlon, self.maxlon, lon_spacing), + labels=[1, 0, 0, 1]) + bm2.drawparallels(np.arange(self.minlat, self.maxlat, lat_spacing), + labels=[1, 0, 0, 1]) bm2.drawcountries() bm2.drawcoastlines() bm2.drawrivers() @@ -744,38 +725,37 @@ def draw_bb_map(self, fig, ax, BB=None, range_rings=None, for nn in range(len(range_rings)): self.plot_range_ring(range_rings[nn], bm=bm2) fig.colorbar(BBmap, ax=ax) - + def draw_beam_terrain_profile(self, fig, ax, ymin=None, ymax=None): '''Draw the Beam height along with terrain and PBB''' if ymin is None: ymin = 0. if ymax is None: ymax = 5000. - - bc, = ax.plot(self.rng_gnd / 1000., self.h / 1000., '-b', - linewidth=3, label='Beam Center') - b3db, = ax.plot(self.rng_gnd / 1000., (self.h + self.a) / 1000., ':b', - linewidth=1.5, label='3 dB Beam width') + + bc, = ax.plot(self.rng_gnd / 1000., self.h / 1000., '-b', + linewidth=3, label='Beam Center') + b3db, = ax.plot(self.rng_gnd / 1000., (self.h + self.a) / 1000., ':b', + linewidth=1.5, label='3 dB Beam width') ax.plot(self.rng_gnd / 1000., (self.h - self.a) / 1000., ':b') - tf = ax.fill_between(self.rng_gnd / 1000., 0., - self.terr[self.Az_plot, :] / 1000., - color='0.75') + tf = ax.fill_between(self.rng_gnd / 1000., 0., + self.terr[self.Az_plot, :] / 1000., color='0.75') ax.set_xlim(0., self.range) ax.set_ylim(ymin / 1000., ymax / 1000.) - ax.set_title(r'dN/dh = %d km$^{-1}$; Elevation angle = %g degrees at %d azimuth'% \ - (self.dNdH, self.E, self.Az_plot), fontdict=TITLEDICT) + ax.set_title(r'dN/dh = %d km$^{-1}$; Elevation angle = %g degrees at %d azimuth' % + (self.dNdH, self.E, self.Az_plot), fontdict=TITLEDICT) ax.set_xlabel('Range (km)') ax.set_ylabel('Height (km)') axb = ax.twinx() - bbf, = axb.plot(self.rng_gnd / 1000., self.CBB[self.Az_plot, :], '-k', - label='BBF') + bbf, = axb.plot(self.rng_gnd / 1000., self.CBB[self.Az_plot, :], '-k', + label='BBF') axb.set_ylabel('Beam-blockage fraction') axb.set_ylim(0., 1.) axb.set_xlim(0., self.range) - + ax.legend((bc, b3db, bbf), ('Beam Center', '3 dB Beam width', 'BBF'), - loc='upper left', fontsize=10) + loc='upper left', fontsize=10) def plot_range_ring(self, range_ring_location_km, bm=None, color='k', ls='-'): @@ -791,19 +771,15 @@ def plot_range_ring(self, range_ring_location_km, bm=None, Axis to plot on. None will use the current axis. """ npts = 100 - bm.tissot(self.rlon, self.rlat, + bm.tissot(self.rlon, self.rlat, np.degrees(range_ring_location_km * 1000. / RE), npts, fill=False, color='black', linestyle='dashed') - -################ -# Save methods # -################ def save_map_data(self, filename=None): """ Save beam blockage and mapped data to NetCDF file - - Parameters:: + + Parameters ---------- filename : str Output name of NetCDF file @@ -812,55 +788,58 @@ def save_map_data(self, filename=None): filename = 'BBmap_data' else: filename = filename - + # Open a NetCDF file to write to nc_fid = nc4.Dataset(filename + '.nc', 'w', format='NETCDF4') nc_fid.description = "PyRadarMet BeamBlockage calculations" - + # Define dimensions xid = nc_fid.createDimension('range', self.nrng) yid = nc_fid.createDimension('azimuth', self.naz) - + # Set global attributes nc_fid.topo_source = 'USGS GTOPO30 DEM files' nc_fid.history = 'Created ' + time.ctime(time.time()) + ' using python netCDF4 module' - + # Create Output variables azid = nc_fid.createVariable('azimuth', np.float32, ('azimuth',)) azid.units = 'degrees' azid.long_name = 'Azimuth' azid[:] = self.phis[:] - + rngid = nc_fid.createVariable('range', np.float32, ('range',)) rngid.units = 'meters' rngid.long_name = 'Range' rngid[:] = self.rng_gnd[:] - - topoid = nc_fid.createVariable('topography', np.float32, ('azimuth','range')) + + topoid = nc_fid.createVariable( + 'topography', np.float32, ('azimuth', 'range')) topoid.units = 'meters' topoid.long_name = 'Terrain Height' topoid[:] = self.terr[:] - - lonid = nc_fid.createVariable('longitude', np.float32, ('azimuth', 'range')) + + lonid = nc_fid.createVariable( + 'longitude', np.float32, ('azimuth', 'range')) lonid.units = 'degrees east' lonid.long_name = 'Longitude' lonid[:] = self.rng_lon[:] - - latid = nc_fid.createVariable('latitude', np.float32, ('azimuth', 'range')) + + latid = nc_fid.createVariable( + 'latitude', np.float32, ('azimuth', 'range')) latid.units = 'degrees north' latid.long_name = 'Latitude' latid[:] = self.rng_lat[:] - - pbbid = nc_fid.createVariable('pbb', np.float32, ('azimuth','range')) + + pbbid = nc_fid.createVariable('pbb', np.float32, ('azimuth', 'range')) pbbid.units = 'unitless' pbbid.long_name = 'Partial Beam Blockage Fraction' pbbid[:] = self.PBB[:] - - cbbid = nc_fid.createVariable('cbb', np.float32, ('azimuth','range')) + + cbbid = nc_fid.createVariable('cbb', np.float32, ('azimuth', 'range')) cbbid.units = 'unitless' cbbid.long_name = 'Cumulative Beam Blockage Fraction' cbbid[:] = self.CBB[:] - + print "Saving " + filename + '.nc' # Close the NetCDF file nc_fid.close() diff --git a/pyradarmet/__init__.py b/pyradarmet/__init__.py index f3866d5..7c24d0f 100644 --- a/pyradarmet/__init__.py +++ b/pyradarmet/__init__.py @@ -1,5 +1,5 @@ """ -PyRadarMet - Python Package that contains a variety of functions for radar meteorology +PyRadarMet - Python Package containing functions for radar meteorology ================================== Top-level package (:mod:`pyradarmet`) ================================== diff --git a/pyradarmet/doppler.py b/pyradarmet/doppler.py index 0251aae..53306f4 100644 --- a/pyradarmet/doppler.py +++ b/pyradarmet/doppler.py @@ -14,7 +14,8 @@ import numpy as np -speed_of_light = 3e8 # Speed of light [m/s] +speed_of_light = 3e8 # Speed of light [m/s] + def freq(lam): """Frequency [Hz] given wavelength. @@ -104,7 +105,8 @@ def Rmax(PRF): def doppler_dilemma(varin, lam): """ The "Doppler dilemma" is the fact that both the Nyquist velocity and - unambiguous maximum range of the radar are based upon the PRF of the system. + unambiguous maximum range of the radar are based upon + the PRF of the system. However, they are inversely proportional, meaning that increasing one requires a decrease in the other. A trade-off inherent in Doppler radar @@ -122,9 +124,10 @@ def doppler_dilemma(varin, lam): """ return (speed_of_light * lam / 8.) / np.asarray(varin) -######################## -## MOBILE PLATFORMS ## -######################## +###################### +# MOBILE PLATFORMS # +###################### + def Vshift(ground_speed, psi): """ @@ -180,4 +183,3 @@ def Vmax_dual(lam, prf1, prf2): Vmax = lam / (4 * ((1. / prf1) - (1. / prf2))) return Vmax - diff --git a/pyradarmet/geometry.py b/pyradarmet/geometry.py index ecd5ab3..0c49427 100644 --- a/pyradarmet/geometry.py +++ b/pyradarmet/geometry.py @@ -15,10 +15,12 @@ import numpy as np -earth_radius = 6371000 # Earth's average radius [m] assuming sphericity -r43 = earth_radius * 4./3. # 4/3 Approximation effective radius for standard atmosphere [m] - # Earth radius taken according to International Union of Geodesy and Geophysics +# Earth's average radius [m] assuming sphericity +earth_radius = 6371000 +# 4/3 Approximation effective radius for standard atmosphere [m] +r43 = earth_radius * 4./3. + def r_effective(dndh=-39e-6): """ @@ -87,7 +89,7 @@ def ray_height(r, elev, h0, reff=r43): """ # Convert earth's radius to km for common dN/dH values and then # multiply by 1000 to return radius in meters - term1 = (np.sqrt(np.asarray(r)**2 +reff**2 + + term1 = (np.sqrt(np.asarray(r)**2 + reff**2 + 2 * np.asarray(r) * reff * np.sin(np.deg2rad(elev)))) h = term1 - reff + h0 return h @@ -216,7 +218,7 @@ def beam_block_frac(Th, Bh, a): # radar beam (Bech et al. (2003), Fig.3) y = Th - Bh - Numer = (y * np.sqrt(a**2 - y**2)) + (a**2 * np.arcsin(y/a)) + (np.pi * a**2 /2.) + Numer = (y * np.sqrt(a**2 - y**2)) + (a**2 * np.arcsin(y/a)) + (np.pi * a**2 / 2.) Denom = np.pi * a**2 diff --git a/pyradarmet/system.py b/pyradarmet/system.py index d0c2295..fcb3b2d 100644 --- a/pyradarmet/system.py +++ b/pyradarmet/system.py @@ -12,8 +12,9 @@ import numpy as np -speed_of_light = 3e8 # Speed of light [m/s] -kBoltz = 1.381e-23 # Boltzmann's constant [ m^2 kg s^-2 K^-1] +speed_of_light = 3e8 # Speed of light [m/s] +kBoltz = 1.381e-23 # Boltzmann's constant [ m^2 kg s^-2 K^-1] + def gain_Pratio(p1, p2): """ @@ -174,7 +175,8 @@ def power_target(power_t, gain, areat, r): def xsec_bscatter_sphere(diam, lam, dielectric=0.93): """ - Backscatter cross-sectional area [m^-2] of a sphere using the Rayleigh approximation. + Backscatter cross-sectional area [m^-2] of a sphere using + the Rayleigh approximation. From Rinehart (1997), Eqn 4.9 and 5.7 @@ -190,9 +192,9 @@ def xsec_bscatter_sphere(diam, lam, dielectric=0.93): Notes ----- The Rayleigh approximation is good when the diamter of a spherical particle - is much smaller than the wavelength of the radar (D/wavelength= 1/16). This - condition leads to the relationship that the area is proportional to the - sixth power of the diameter. + is much smaller than the wavelength of the radar (D/wavelength= 1/16). + This condition leads to the relationship that the area is proportional + to the sixth power of the diameter. The default is for a dielectric factor value for water. This can be changed by the user, e.g. K=0.208 for particle sizes of equivalent melted @@ -203,7 +205,8 @@ def xsec_bscatter_sphere(diam, lam, dielectric=0.93): def norm_xsec_bscatter_sphere(diam, lam, dielectric=0.93): """ - Normalized Backscatter cross-sectional area [m^2] of a sphere using the Rayleigh approximation. + Normalized Backscatter cross-sectional area [m^2] of a sphere using + the Rayleigh approximation. From Rinehart (1997), Eqn 4.9 and 5.7 and Battan Ch. 4.5 @@ -219,9 +222,9 @@ def norm_xsec_bscatter_sphere(diam, lam, dielectric=0.93): Notes ----- The Rayleigh approximation is good when the diamter of a spherical particle - is much smaller than the wavelength of the radar (D/wavelength= 1/16). This - condition leads to the relationship that the area is proportional to the - sixth power of the diameter. + is much smaller than the wavelength of the radar (D/wavelength= 1/16). + This condition leads to the relationship that the area is proportional + to the sixth power of the diameter. The default is for a dielectric factor value for water. This can be changed by the user, e.g. K=0.208 for particle sizes of equivalent melted @@ -230,7 +233,7 @@ def norm_xsec_bscatter_sphere(diam, lam, dielectric=0.93): # Calculate the cross-sectional backscatter area sig = xsec_bscatter_sphere(np.asarray(diam), lam, dielectric) - return sig/ (np.pi * (np.asarray(diam)/2.)**2) + return sig / (np.pi * (np.asarray(diam)/2.)**2) def size_param(diam, lam): @@ -248,8 +251,9 @@ def size_param(diam, lam): Notes ----- - The size paramter can be used along with the backscattering cross-section to - distinguish ice and water dielectric characteristics. For example: + The size paramter can be used along with the backscattering + cross-section todistinguish ice and water dielectric characteristics. + For example: Alpha < 2 the backscattering cross-section of ice is smaller than water, Alpha > 2 the opposite is true due to the fact that absorption in water exceeds that in ice. @@ -259,7 +263,8 @@ def size_param(diam, lam): def power_return_target(power_t, gain, lam, sig, r): """ - Power [W] returned y target located at the center of the antenna beam pattern. + Power [W] returned y target located at the center + of the antenna beam pattern. From Rinehart (1997), Eqn 4.7 @@ -304,9 +309,9 @@ def thermal_noise(bandwidth, Units, noise_temp=290.): # Calculate the noise, convert if requested noise = kBoltz * noise_temp * np.asarray(bandwidth) - if Units.upper()=='W': + if Units.upper() == 'W': noiset = noise - elif Units.upper()=='DBM': + elif Units.upper() == 'DBM': noiset = 10. * np.log10(noise/10**-3) else: print("Units must be in 'W' or 'dBm'") diff --git a/pyradarmet/variables.py b/pyradarmet/variables.py index 63579d7..f0748e1 100644 --- a/pyradarmet/variables.py +++ b/pyradarmet/variables.py @@ -48,11 +48,12 @@ def reflectivity(power_t, gain, pulse_width, wavelength, bw_h, bw_v, ----- This routine calls the radar_constant function. The default is for a dielectric factor value for water. This can be - changed by the user, e.g. K=0.208 for particle sizes of equivalent melted - diameters or K=0.176 for particle sizes of equivalent ice spheres. + changed by the user, e.g. K=0.208 for particle sizes of equivalent melted + diameters or K=0.176 for particle sizes of equivalent ice spheres. """ # Call the radar constant function - C1 = radar_constant(power_t, gain, pulse_width, wavelength, bw_h, bw_v, aloss, rloss) + C1 = radar_constant(power_t, gain, pulse_width, wavelength, + bw_h, bw_v, aloss, rloss) return power_return * np.asarray(r)**2 / (C1 * dielectric**2) @@ -191,7 +192,7 @@ def zdp(z_h, z_v): zdp = np.full_like(zh, np.nan) good = np.where(zh > zv) - zdp[good] = 10.* np.log10(zh[good] - zv[good]) + zdp[good] = 10. * np.log10(zh[good] - zv[good]) return zdp @@ -217,7 +218,8 @@ def hdr(dbz_h, zdr): Considerations for this equation (see paper for more details): 1) Standar error of disdrometer data allowed for - 2) Drop oscillation accounted for based on 50% model of Seliga et al (1984) + 2) Drop oscillation accounted for based on 50% model + of Seliga et al (1984) 3) Lower (27) bound chose to provide constant Zh ref level 4) Upper cutoff of 60 (may be too low) diff --git a/pyradarmet/zdrcal.py b/pyradarmet/zdrcal.py index 4837c5e..cb6e33c 100644 --- a/pyradarmet/zdrcal.py +++ b/pyradarmet/zdrcal.py @@ -1,13 +1,15 @@ """ pyradarmet.zdrcal -========================= +================= Calculate the differential reflectivity (Zdr) offset from polarimetric radars. -Adapted by Nick Guy from original Fortran code written by David Jorgensen, - implemented more robust processing to reduce noise in data. Thanks to Joseph - Hardin regarding processing and posting a notebook online that got this started. - http://www.josephhardinee.com/blog/?p=35 +Adapted by Nick Guy from original Fortran code written by David Jorgensen. +Implemented a more robust processing to reduce noise in data. +Thanks to Joseph Hardin regarding processing and posting a notebook online +that got this started. + +http://www.josephhardinee.com/blog/?p=35 .. autosummary:: @@ -15,27 +17,26 @@ calculate_zdr_offset - """ import numpy as np import matplotlib.pylab as plt import pyart -#from .conversion import dBZ2Z, Z2dBZ -#====================================================================== -def calculate_zdr_offset(radar, debug=False, - remove_first_n_gates=None, sample_min=100, - rhv_min=None, sig_min=0.5,Htmax=None, dbz_min=None, - zdr_field=None, refl_field=None, phidp_field=None, - rhv_field=None, kdp_field=None): + + +def calculate_zdr_offset(radar, debug=False, + remove_first_n_gates=None, sample_min=100, + rhv_min=None, sig_min=0.5, Htmax=None, dbz_min=None, + zdr_field=None, refl_field=None, phidp_field=None, + rhv_field=None, kdp_field=None): """ Differential reflectivity calibration from 'birdbath' scans. Parameters ---------- radar : Radar - Radar object (from PyArt) to use for attenuation calculations. Must have - copol_coeff, norm_coherent_power, proc_dp_phase_shift, + Radar object (from PyArt) to use for attenuation calculations. + Must have copol_coeff, norm_coherent_power, proc_dp_phase_shift, reflectivity_horizontal fields. debug : bool True to print debugging information, False supressed this printing. @@ -55,7 +56,8 @@ def calculate_zdr_offset(radar, debug=False, rhv_min : float Minimum copol_coeff value to consider valid. sig_min : float - Minimum standard deviation to consider valid and include in calculations. + Minimum standard deviation to consider valid and include + in calculations. dbz_min : float Minimum reflectivity to consider valid and include in calculations. Htmax : float @@ -69,19 +71,18 @@ def calculate_zdr_offset(radar, debug=False, Usage ---------- - """ - # Pull out some information to send back in zdr_bias dictionary - Rlat = radar.latitude['data'] # Radar latitude - Rlon = radar.longitude['data'] # Radar longitude - RAlt = radar.altitude['data'] # Radar altitude - RadName = radar.metadata['instrument_name'] # Radar name - SitName = radar.metadata['source'] # Site/project name - GenName = radar.metadata['institution'] # Facilitiy name - time = radar.time # Time at middle of each ray - range = radar.range # Range along ray - azimuth = radar.azimuth # Azimuth of scan - + # Pull out some information to send back in zdr_bias dictionary + Rlat = radar.latitude['data'] # Radar latitude + Rlon = radar.longitude['data'] # Radar longitude + RAlt = radar.altitude['data'] # Radar altitude + RadName = radar.metadata['instrument_name'] # Radar name + SitName = radar.metadata['source'] # Site/project name + GenName = radar.metadata['institution'] # Facilitiy name + time = radar.time # Time at middle of each ray + rng = radar.range # Range along ray + azimuth = radar.azimuth # Azimuth of scan + # Create 2D arrays for masking data later # Ht2D, az2D = np.meshgrid(range['data']/1000.,azimuth['data']) # print radar.fields.keys() @@ -89,15 +90,15 @@ def calculate_zdr_offset(radar, debug=False, if zdr_field is None: zdr_field = pyart.config.get_field_name('differential_reflectivity') if zdr_field not in radar.fields: - zdr_field = pyart.config.get_field_name('ZDR') + zdr_field = pyart.config.get_field_name('ZDR') if refl_field is None: refl_field = pyart.config.get_field_name('reflectivity') if refl_field not in radar.fields: - refl_field = pyart.config.get_field_name('DBZ') + refl_field = pyart.config.get_field_name('DBZ') if rhv_field is None: rhv_field = pyart.config.get_field_name('cross_correlation_ratio') if rhv_field not in radar.fields: - rhv_field = pyart.config.get_field_name('RHOHV') + rhv_field = pyart.config.get_field_name('RHOHV') if phidp_field is None: # use corrrected_differential_phae or unfolded_differential_phase # fields if they are available, if not use differential_phase field @@ -117,45 +118,44 @@ def calculate_zdr_offset(radar, debug=False, Rhohv = radar.fields[rhv_field]['data'] PhiDP = radar.fields[phidp_field]['data'] KDP = radar.fields[kdp_field]['data'] - + # Extract parameters from radar # nsweeps = int(radar.nsweeps) - nGates = int(radar.ngates) # Number of gates - nrays = int(radar.nrays) # Number of rays + nGates = int(radar.ngates) # Number of gates + nrays = int(radar.nrays) # Number of rays # Apply mask across variables for missing data ZDR, dBZ, Rhv, PhiDP, KDP = mask_variables(ZDR, dBZ, Rhohv, PhiDP, KDP) - # Find corresponding axes for range (height) and azimuth + # Find corresponding axes for range (height) and azimuth axAz, axHt = get_ray_range_dims(ZDR, nrays) - + # Apply mask for chosen gates nearest to radar and above height threshold - ZDR, dBZ, Rhv, PhiDP, KDP = mask_gates(ZDR, dBZ, Rhv, PhiDP,KDP, - range['data']/1000., azimuth['data'], axAz, - Htmax=Htmax, remove_first_n_gates=remove_first_n_gates) - + ZDR, dBZ, Rhv, PhiDP, KDP = mask_gates( + ZDR, dBZ, Rhv, PhiDP, KDP, rng['data']/1000., azimuth['data'], axAz, + Htmax=Htmax, remove_first_n_gates=remove_first_n_gates) + # Apply mask for reflectivity below a threshold (proxy for low SNR) ZDR, dBZ, Rhv, PhiDP, KDP = mask_refl(ZDR, dBZ, Rhv, PhiDP, KDP, mindBZ=dbz_min) - # Find individual ray properties (along each ray at azimuthal step) DR_RaySum, DR_RayAvg, DR_RayStd, NGood_Ray = calc_stats(ZDR, axAz) - - # Find individual range properties (along a constant distance [height] from radar) + + # Find individual range properties at a constant distance [height] from radar DR_HtSum, DR_HtAvg, DR_HtStd, NGood_Ht = calc_stats(ZDR, axHt) RH_HtSum, RH_HtAvg, RH_HtStd, NGood_RH_Ht = calc_stats(Rhv, axHt) dBZ_HtSum, dBZ_HtAvg, dBZ_HtStd, NGood_dBZ_Ht = calc_stats(dBZ, axHt) KDP_HtSum, KDP_HtAvg, KDP_HtStd, NGood_KDP_Ht = calc_stats(KDP, axHt) - # Apply mask for data with range averaged Std Dev and azimuthal averaged Std Dev - # greater than threshold - ZDR, dBZ, Rhv, PhiDP, KDP = mask_sigma(ZDR, dBZ, Rhv, PhiDP, KDP, DR_HtStd, DR_RayStd, - axAz, minSig=sig_min) - - # Find volume properites + # Apply mask for data with range averaged Std Dev and azimuthal + # averaged Std Deviations greater than threshold + ZDR, dBZ, Rhv, PhiDP, KDP = mask_sigma( + ZDR, dBZ, Rhv, PhiDP, KDP, DR_HtStd, DR_RayStd, axAz, minSig=sig_min) + + # Find volume properites DR_Avg = np.ma.mean(ZDR) DR_Std = np.ma.std(ZDR, ddof=1) NGood = np.ma.count(ZDR) - + if (NGood > sample_min): # Find statistics for filtered (masked) variables DR_HtSum_filt, DR_HtAvg_filt, DR_HtStd_filt, NGood_Ht_filt = calc_stats(ZDR, axHt) @@ -163,42 +163,42 @@ def calculate_zdr_offset(radar, debug=False, dBZ_HtSum_filt, dBZ_HtAvg_filt, dBZ_HtStd_filt, NGood_dBZ_Ht_filt = calc_stats(dBZ, axHt) KDP_HtSum_filt, KDP_HtAvg_filt, KDP_HtStd_filt, NGood_KDP_Ht_filt = calc_stats(KDP, axHt) PhiDP_HtSum_filt, PhiDP_HtAvg_filt, PhiDP_HtStd_filt, NGood_PhiDP_Ht_filt = calc_stats(PhiDP, axHt) - + # Create an output dictionary - zdr_bias = create_dict(DR_Avg, DR_Std, NGood, - DR_RaySum, DR_RayAvg, DR_RayStd, NGood_Ray, - DR_HtAvg, DR_HtStd, NGood_Ht, - RH_HtAvg, dBZ_HtAvg, KDP_HtAvg, - ZDR, DR_HtAvg_filt, DR_HtStd_filt, - RH_HtAvg_filt, dBZ_HtAvg_filt, KDP_HtAvg_filt, PhiDP_HtAvg_filt, - RadName, GenName, Rlat, Rlon, RAlt, time, range, azimuth) - + zdr_bias = create_dict( + DR_Avg, DR_Std, NGood, DR_RaySum, DR_RayAvg, DR_RayStd, NGood_Ray, + DR_HtAvg, DR_HtStd, NGood_Ht, RH_HtAvg, dBZ_HtAvg, KDP_HtAvg, + ZDR, DR_HtAvg_filt, DR_HtStd_filt, RH_HtAvg_filt, dBZ_HtAvg_filt, + KDP_HtAvg_filt, PhiDP_HtAvg_filt, RadName, GenName, Rlat, Rlon, + RAlt, time, rng, azimuth) + else: zdr_bias = None # Send back the dictionary containing data. return zdr_bias -#====================================================================== + + def mask_variables(ZDR, dBZ, RhoHV, PhiDP, KDP, rhv_min=0.8): # Combine the masks for all variables to ensure no missing data points maskZDR = np.ma.getmask(ZDR) maskZ = np.ma.getmask(dBZ) maskRhv = np.ma.getmask(RhoHV) maskPdp = np.ma.getmask(PhiDP) - is_cor = RhoHV < rhv_min # Mask values where Correl Coeff < threshold - - mskt1 = np.logical_and(maskZ, maskZDR) # Combine refl and ZDR masks - mskt2 = np.logical_and(mskt1, maskRhv) # Combine with missing Rho HV mask - mskt3 = np.logical_and(mskt2, is_cor) # Combine with min threshold Rho HV - mask = np.logical_and(mskt3, maskPdp) # Combine with missing Phidp mask - + is_cor = RhoHV < rhv_min # Mask values where Correl Coeff < threshold + + mskt1 = np.logical_and(maskZ, maskZDR) # Combine refl and ZDR masks + mskt2 = np.logical_and(mskt1, maskRhv) # Combine with missing Rho HV mask + mskt3 = np.logical_and(mskt2, is_cor) # Combine with min threshold Rho HV + mask = np.logical_and(mskt3, maskPdp) # Combine with missing Phidp mask + ZDR = np.ma.masked_where(mask, ZDR) dBZ = np.ma.masked_where(mask, dBZ) Rhv = np.ma.masked_where(mask, RhoHV) PhiDP = np.ma.masked_where(mask, PhiDP) KDP = np.ma.masked_where(mask, KDP) - return ZDR, dBZ, Rhv, PhiDP, KDP -#====================================================================== + + def get_ray_range_dims(ZDR, nrays): if nrays == ZDR.shape[0]: axAz = 1 @@ -210,22 +210,23 @@ def get_ray_range_dims(ZDR, nrays): print "Makes sure you are using polar coordinate data!" return return axAz, axHt -#====================================================================== -def mask_gates(ZDR, dBZ, Rhv, PhiDP, KDP, + + +def mask_gates(ZDR, dBZ, Rhv, PhiDP, KDP, Ht, Az, axAz, Htmax=10., remove_first_n_gates=None): # Create 2D arrays for masking data later Ht2D, az2D = np.meshgrid(Ht, Az) - + # Check to see what gate to start with user can set or default to 0 if remove_first_n_gates is None: - GateBeg=1 + GateBeg = 1 elif remove_first_n_gates == 0: - GateBeg=1 + GateBeg = 1 else: - GateBeg=remove_first_n_gates # because of 0 based indexing + GateBeg = remove_first_n_gates # because of 0 based indexing # Check which axis to perform calculations over - # Mask out the lower gates based upon remove_first_n_gates keyword + # Mask out the lower gates based upon remove_first_n_gates keyword if axAz == 1: # Mask lower gates ZDR[:, 0:GateBeg-1] = np.ma.masked @@ -233,7 +234,7 @@ def mask_gates(ZDR, dBZ, Rhv, PhiDP, KDP, Rhv[:, 0:GateBeg-1] = np.ma.masked PhiDP[:, 0:GateBeg-1] = np.ma.masked KDP[:, 0:GateBeg-1] = np.ma.masked - + # Mask upper gates ZDR = np.ma.masked_where(Ht2D > Htmax, ZDR) dBZ = np.ma.masked_where(Ht2D > Htmax, dBZ) @@ -246,23 +247,22 @@ def mask_gates(ZDR, dBZ, Rhv, PhiDP, KDP, Rhv[0:GateBeg-1, :] = np.ma.masked PhiDP[0:GateBeg-1, :] = np.ma.masked KDP[0:GateBeg-1, :] = np.ma.masked - + # Mask upper gates ZDR = np.ma.masked_where(Ht2D.transpose() > Htmax, ZDR) dBZ = np.ma.masked_where(Ht2D.transpose() > Htmax, dBZ) Rhv = np.ma.masked_where(Ht2D.transpose() > Htmax, Rhv) PhiDP = np.ma.masked_where(Ht2D.transpose() > Htmax, PhiDP) KDP = np.ma.masked_where(Ht2D.transpose() > Htmax, KDP) - return ZDR, dBZ, Rhv, PhiDP, KDP -#====================================================================== -def mask_sigma(ZDR, dBZ, Rhv, PhiDP, KDP, HtStd, AzStd, axAz, minSig=0.5): + +def mask_sigma(ZDR, dBZ, Rhv, PhiDP, KDP, HtStd, AzStd, axAz, minSig=0.5): # Create 2D arrays for masking data later HtStd2D, RayStd2D = np.meshgrid(HtStd, AzStd) # Check which axis to perform calculations over - # Mask out the lower gates based upon remove_first_n_gates keyword + # Mask out the lower gates based upon remove_first_n_gates keyword if axAz == 1: # Mask along the range (height) ZDR = np.ma.masked_where(HtStd2D > minSig, ZDR) @@ -270,7 +270,7 @@ def mask_sigma(ZDR, dBZ, Rhv, PhiDP, KDP, HtStd, AzStd, axAz, minSig=0.5): Rhv = np.ma.masked_where(HtStd2D > minSig, Rhv) PhiDP = np.ma.masked_where(HtStd2D > minSig, PhiDP) KDP = np.ma.masked_where(HtStd2D > minSig, KDP) - + # Mask along the Azimuths ZDR = np.ma.masked_where(RayStd2D > minSig, ZDR) dBZ = np.ma.masked_where(RayStd2D > minSig, dBZ) @@ -284,16 +284,16 @@ def mask_sigma(ZDR, dBZ, Rhv, PhiDP, KDP, HtStd, AzStd, axAz, minSig=0.5): Rhv = np.ma.masked_where(HtStd2D.transpose() > minSig, Rhv) PhiDP = np.ma.masked_where(HtStd2D.transpose() > minSig, PhiDP) KDP = np.ma.masked_where(HtStd2D.transpose() > minSig, KDP) - + # Mask along the Azimuths ZDR = np.ma.masked_where(RayStd2D.transpose() > minSig, ZDR) dBZ = np.ma.masked_where(RayStd2D.transpose() > minSig, dBZ) Rhv = np.ma.masked_where(RayStd2D.transpose() > minSig, Rhv) PhiDP = np.ma.masked_where(RayStd2D.transpose() > minSig, PhiDP) KDP = np.ma.masked_where(RayStd2D.transpose() > minSig, KDP) - return ZDR, dBZ, Rhv, PhiDP, KDP -#====================================================================== + + def mask_refl(ZDR, dBZ, Rhv, PhiDP, KDP, mindBZ=-40.): # Mask out the points with reflectivity < mindBZ ZDR = np.ma.masked_where(dBZ < mindBZ, ZDR) @@ -301,71 +301,72 @@ def mask_refl(ZDR, dBZ, Rhv, PhiDP, KDP, mindBZ=-40.): Rhv = np.ma.masked_where(dBZ < mindBZ, Rhv) PhiDP = np.ma.masked_where(dBZ < mindBZ, PhiDP) KDP = np.ma.masked_where(dBZ < mindBZ, KDP) - return ZDR, dBZ, Rhv, PhiDP, KDP -#====================================================================== + + def calc_stats(Var, Ax): Sum = np.ma.cumsum(Var, axis=Ax) Avg = np.ma.mean(Var, axis=Ax) Std = np.ma.std(Var, axis=Ax) nValid = np.ma.count(Var, axis=Ax) - return Sum, Avg, Std, nValid -#====================================================================== -def create_dict(DR_Avg, DR_Std, NGood, + + +def create_dict(DR_Avg, DR_Std, NGood, DR_RaySum, DR_RayAvg, DR_RayStd, NGood_Ray, DR_HtAvg, DR_HtStd, NGood_Ht, RH_HtAvg, dBZ_HtAvg, KDP_HtAvg, ZDR, DR_HtAvg_filt, DR_HtStd_filt, - RH_HtAvg_filt, dBZ_HtAvg_filt, KDP_HtAvg_filt, PhiDP_HtAvg_filt, - RadName, GenName, Rlat, Rlon, RAlt, time, range, azimuth): - - zdr_bias = {'volume_average': DR_Avg, + RH_HtAvg_filt, dBZ_HtAvg_filt, KDP_HtAvg_filt, + PhiDP_HtAvg_filt, RadName, GenName, Rlat, Rlon, RAlt, + time, rng, azimuth): + + zdr_bias = {'volume_average': DR_Avg, 'volume_standard_deviation': DR_Std, 'volume_number_good': NGood, - 'ray_cumulative_sum': DR_RaySum, + 'ray_cumulative_sum': DR_RaySum, 'ray_average': DR_RayAvg, - 'ray_standard_deviation': DR_RayStd, + 'ray_standard_deviation': DR_RayStd, 'ray_number_good': NGood_Ray, 'height_average': DR_HtAvg, - 'height_standard_deviation': DR_HtStd, - 'height_number_good': NGood_Ht, - 'rhv_height_average': RH_HtAvg, + 'height_standard_deviation': DR_HtStd, + 'height_number_good': NGood_Ht, + 'rhv_height_average': RH_HtAvg, 'dbz_height_average': dBZ_HtAvg, 'kdp_height_average': KDP_HtAvg, - 'filtered_zdr': ZDR, + 'filtered_zdr': ZDR, 'height_average_filt': DR_HtAvg_filt, 'height_standard_deviation_filt': DR_HtStd_filt, - 'rhv_height_average_filt': RH_HtAvg_filt, + 'rhv_height_average_filt': RH_HtAvg_filt, 'dbz_height_average_filt': dBZ_HtAvg_filt, 'kdp_height_average_filt': KDP_HtAvg_filt, 'phidp_height_average_filt': PhiDP_HtAvg_filt, - 'radar_name': RadName, + 'radar_name': RadName, 'facility_name': GenName, - 'radar_latitude':Rlat, - 'radar_longitude':Rlon, - 'radar_altitude':RAlt, - 'time': time, - 'range': range, - 'azimuth':azimuth} - + 'radar_latitude': Rlat, + 'radar_longitude': Rlon, + 'radar_altitude': RAlt, + 'time': time, + 'range': rng, + 'azimuth': azimuth} return zdr_bias -#====================================================================== + + def plot_zdr_rhv(X1, X2, Y, xlims1=None, xlims2=None, ymax=None, - label=True, labelTx=' ', labelpos=None): + label=True, labelTx=' ', labelpos=None): """Create a vertical plot of differential reflectivity and copolar correlation coefficient. """ fig, ax = plt.subplots() axes = [ax, ax.twiny()] - + if xlims1 is None: - xlims1 = (-1.5, 1.5) + xlims1 = (-1.5, 1.5) if xlims2 is None: - xlims2 = (.8, 1.) + xlims2 = (.8, 1.) if ymax is None: - ymax = 10. - + ymax = 10. + axes[0].plot(X1, Y, label='Differential Reflectivity', color='b') axes[1].plot(X2, Y, label='CoPolar Correlation Coefficient', color='r') axes[0].set_xlim(xlims1) @@ -382,12 +383,13 @@ def plot_zdr_rhv(X1, X2, Y, xlims1=None, xlims2=None, ymax=None, axes[0].xaxis.grid(color='b', ls=':', lw=1) axes[1].xaxis.grid(color='r', ls=':', lw=1) axes[0].yaxis.grid(color='k', ls=':', lw=.5) - + handles1, labels1 = axes[0].get_legend_handles_labels() handles2, labels2 = axes[1].get_legend_handles_labels() plt.legend([handles1[0], handles2[0]], - ['Differential Reflectivity','CoPolar Correlation Coefficient'], - loc=6, fontsize =11) + ['Differential Reflectivity', + 'CoPolar Correlation Coefficient'], + loc=6, fontsize=11) # Add label to lower left in plot unless location specified if label: if labelpos is None: @@ -395,46 +397,49 @@ def plot_zdr_rhv(X1, X2, Y, xlims1=None, xlims2=None, ymax=None, plt.figtext(labelpos[0], labelpos[1], labelTx) return fig, ax, axes -#====================================================================== + + def oplot_zdr_rhv(X1, X2, Y, ax=None, alf=1.): """Add additional differential reflectivity and copolar correlation coefficient lines to existing graphic.""" ax[0].plot(X1, Y, color='b', alpha=alf) ax[1].plot(X2, Y, color='r', alpha=alf) -#====================================================================== + + def plot_add_stats(instance=None, avg=None, sd=None, points=None, good_samples=None, tot_samples=None, labelpos=None): """Add text to the figure that contains statistics.""" if labelpos is None: - labelpos=(.15, .7) + labelpos = (.15, .7) if instance == 'single': plt.figtext(labelpos[0], labelpos[1], - "Avg Zdr = %g\nStd Dev = %g\nPoints = %g\n"%(avg, sd, points)) + "Avg Zdr = %g\nStd Dev = %g\nPoints = %g\n" % (avg, sd, points)) elif instance == 'multi': plt.figtext(labelpos[0], labelpos[1], - "Avg Zdr = %g\nStd Dev = %g\nSamples = %g used of %g total\nPoints = %g\n"% - (avg, sd, good_samples, tot_samples, points)) + "Avg Zdr = %g\nStd Dev = %g\nSamples = %g used of %g total\nPoints = %g\n" % + (avg, sd, good_samples, tot_samples, points)) else: - print "Need to supply either single or multiple run statistics!!" - -#====================================================================== -def plot_3panel(X1, X2, X3, Y, - xlims1=None, xlims2=None, xlims3=None, ymax=None, - X1lab='dBZ', X2lab='ZDR', X3lab=r'Kdp (dB km$^{-1}$)', ylab=None, - label=True, labelTx=' ', labelpos=None): - """Create a 3-panel vertical plot of reflectivity, differential reflectivity, - and specific differential phase. + print("Need to supply either single or multiple run statistics!!") + + +def plot_3panel(X1, X2, X3, Y, + xlims1=None, xlims2=None, xlims3=None, ymax=None, + X1lab='dBZ', X2lab='ZDR', X3lab=r'Kdp (dB km$^{-1}$)', + ylab=None, label=True, labelTx=' ', labelpos=None): + """ + Create a 3-panel vertical plot of reflectivity, differential reflectivity, + and specific differential phase. """ fig, (ax1, ax2, vax3) = plt.subplots(1, 3, sharey=True) - + if xlims1 is None: - xlims1 = (-20, 30.) + xlims1 = (-20, 30.) if xlims2 is None: - xlims2 = (-1.5, 1.5) + xlims2 = (-1.5, 1.5) if xlims3 is None: - xlims3 = (.2, .4) + xlims3 = (.2, .4) if ymax is None: - ymax = 10. + ymax = 10. ax1.plot(X1, Y, color='k') ax2.plot(X2, Y, color='k') ax3.plot(X3, Y, color='k') @@ -444,15 +449,15 @@ def plot_3panel(X1, X2, X3, Y, ax1.set_xlabel(X1lab) ax2.set_xlabel(X2lab) ax3.set_xlabel(X3lab) - + # Set y axis characteristics ax1.set_ylim(0., ymax) if ylab is None: ax1.set_ylabel('Altitude (km)') else: ax1.set_ylabel(ylab) - ax1.vlines(0,0,ymax) - ax2.vlines(0,0,ymax) + ax1.vlines(0, 0, ymax) + ax2.vlines(0, 0, ymax) ax1.xaxis.grid(color='k', ls=':', lw=1) ax2.xaxis.grid(color='k', ls=':', lw=1) ax3.xaxis.grid(color='k', ls=':', lw=1) @@ -462,15 +467,15 @@ def plot_3panel(X1, X2, X3, Y, # Add label to above plot unless location specified if label: if labelpos is None: - labelpos = (.15,.91) + labelpos = (.15, .91) plt.figtext(labelpos[0], labelpos[1], labelTx) return fig, ax1, ax2, ax3 -#====================================================================== + + def oplot_3panel(X1, X2, X3, Y, ax1=None, ax2=None, ax3=None, alf=1.): """Add additional differential reflectivity and copolar correlation coefficient lines to existing graphic.""" ax1.plot(X1, Y, color='k', alpha=alf) ax2.plot(X2, Y, color='k', alpha=alf) ax3.plot(X3, Y, color='k', alpha=alf) -#====================================================================== \ No newline at end of file diff --git a/scripts/cal_zdr b/scripts/cal_zdr index df166a7..c4138ff 100755 --- a/scripts/cal_zdr +++ b/scripts/cal_zdr @@ -1,7 +1,5 @@ #! /usr/bin/env python -#====================================================================== -# Import needed libraries import pyart import argparse import numpy as np @@ -10,23 +8,24 @@ import matplotlib.pylab as plt from glob import glob from pyradarmet.zdrcal import calculate_zdr_offset as zdrcal -from pyradarmet.zdrcal import plot_zdr_rhv,plot_add_stats,oplot_zdr_rhv -#====================================================================== +from pyradarmet.zdrcal import (plot_zdr_rhv, plot_add_stats, oplot_zdr_rhv) + + # The user can modify these next variables to control how the program processes data -MINSIG = 0.5 # Maximum standard deviation to accept in analysis -MINSAMPLE = 100 # Minimum number of samples used to perform calculations -MINDBZ = 5. # Minimum reflectivity to accept in analysis (proxy for SNR) -GATESKIP = 7 # The number of gates at the beginning or array (closest to radar) to cut out -MINRHV = 0.85 # Minimum correlation coefficient to accept in analysis -MAXHT = 5. # Maximum height to accept in analysis -ZH_LIMS = [-20., 30.] # Range of Reflectivity to plot -RH_LIMS = [0.8, 1.0] # Range of Copolar Correlation Coefficient to plot -DR_LIMS = [-1.5, 1.5] # Range of Differential reflectivity to plot -KD_LIMS = [0.2, 0.4] # Range of Specific Differential Phase to plot +MINSIG = 0.5 # Maximum standard deviation to accept in analysis +MINSAMPLE = 100 # Minimum number of samples used to perform calculations +MINDBZ = 5. # Minimum reflectivity to accept in analysis (proxy for SNR) +GATESKIP = 7 # The number of gates at the beginning or array (closest to radar) to cut out +MINRHV = 0.85 # Minimum correlation coefficient to accept in analysis +MAXHT = 5. # Maximum height to accept in analysis +ZH_LIMS = [-20., 30.] # Range of Reflectivity to plot +RH_LIMS = [0.8, 1.0] # Range of Copolar Correlation Coefficient to plot +DR_LIMS = [-1.5, 1.5] # Range of Differential reflectivity to plot +KD_LIMS = [0.2, 0.4] # Range of Specific Differential Phase to plot pOut = "zdr_bias" pType = "png" -#====================================================================== + if __name__ == '__main__': @@ -58,27 +57,26 @@ if __name__ == '__main__': parser.add_argument('-v', '--version', action='version', version='Py-ART version %s' % (pyart.__version__)) - + parser.add_argument('-p', '--plot', action='store_true', help='create an output plot file') # parser.add_argument('-t', '--text', action='store_true', # help='create and output text file') - + args = parser.parse_args() - # ====================================================================== - + # Search for the file(s) flist = glob(args.searchstring) - + multiStat = False if len(flist) > 1: multiStat = True # Create an array to hold data, only really used if multiple files - Stats = np.ma.empty([len(flist),3]) + Stats = np.ma.empty([len(flist), 3]) index = 0 indexGood = 0 - OutTx = open('zdr_offset.out','w') + OutTx = open('zdr_offset.out', 'w') for file in flist: fName, fileExt = os.path.splitext(file) # read in the file @@ -98,83 +96,83 @@ if __name__ == '__main__': radar = pyart.io.read(file) # Call the bias calculation routine - zdr_bias = zdrcal(radar, remove_first_n_gates=GATESKIP, sample_min=MINSAMPLE,\ - rhv_min=MINRHV, sig_min=MINSIG, dbz_min=MINDBZ, Htmax=MAXHT) - + zdr_bias = zdrcal(radar, remove_first_n_gates=GATESKIP, sample_min=MINSAMPLE, + rhv_min=MINRHV, sig_min=MINSIG, dbz_min=MINDBZ, + Htmax=MAXHT) + # If minimun number of samples is met if zdr_bias is not None: - - txVar = ("Avg Zdr: "+str(zdr_bias['volume_average'])+\ - " Std Dev: "+str(zdr_bias['volume_standard_deviation'])+\ - " Good Pts: "+str(zdr_bias['volume_number_good'])+\ - " File: "+os.path.basename(fName)) - + txVar = ("Avg Zdr: " + str(zdr_bias['volume_average']) + + " Std Dev: " + str(zdr_bias['volume_standard_deviation']) + + " Good Pts: " + str(zdr_bias['volume_number_good']) + + " File: " + os.path.basename(fName)) + if zdr_bias['volume_standard_deviation'] > MINSIG: - Stats[index,0] = np.ma.masked - Stats[index,1] = np.ma.masked - Stats[index,2] = np.ma.masked + Stats[index, 0] = np.ma.masked + Stats[index, 1] = np.ma.masked + Stats[index, 2] = np.ma.masked else: - Stats[index,0] = zdr_bias['volume_average'] - Stats[index,1] = zdr_bias['volume_standard_deviation'] - Stats[index,2] = zdr_bias['volume_number_good'] - + Stats[index, 0] = zdr_bias['volume_average'] + Stats[index, 1] = zdr_bias['volume_standard_deviation'] + Stats[index, 2] = zdr_bias['volume_number_good'] + # Create plot if desired if args.plot: if index == 0: - fig, ax, axes = plot_zdr_rhv(zdr_bias['height_average'],\ - zdr_bias['rhv_height_average'],\ - zdr_bias['range']['data']/1000.,\ - ymax=MAXHT,\ - xlims1=DR_LIMS, xlims2=RH_LIMS,\ - label=True, labelTx=zdr_bias['radar_name']) + fig, ax, axes = plot_zdr_rhv( + zdr_bias['height_average'], + zdr_bias['rhv_height_average'], + zdr_bias['range']['data']/1000., + ymax=MAXHT, xlims1=DR_LIMS, xlims2=RH_LIMS, + label=True, labelTx=zdr_bias['radar_name']) else: -## if zdr_bias['volume_standard_deviation'] > MINSIG: +# if zdr_bias['volume_standard_deviation'] > MINSIG: if zdr_bias['height_standard_deviation'].mean() > MINSIG: alff = 0.1 else: alff = 1. - oplot_zdr_rhv(zdr_bias['height_average'],\ - zdr_bias['rhv_height_average'],\ - zdr_bias['range']['data']/1000., \ + oplot_zdr_rhv(zdr_bias['height_average'], + zdr_bias['rhv_height_average'], + zdr_bias['range']['data']/1000., ax=axes, alf=alff) # Increment indices index += 1 -## if zdr_bias['volume_standard_deviation'] <= MINSIG: +# if zdr_bias['volume_standard_deviation'] <= MINSIG: if zdr_bias['height_standard_deviation'].mean() <= MINSIG: indexGood += 1 - + RadName = zdr_bias['radar_name'] - + del zdr_bias # If zdr_bias returns None, then too few samples were uncovered. else: txVar = ("Not enough samples in volume") RadName = "No_Data" - - print txVar - OutTx.write(txVar+"\n") - - OutTx.write("ALL FILES WITH STD DEV < "+str(MINSIG)+" YIELD ("+str(indexGood)+" files)::\n") - txVar = ("Avg Zdr: "+str(Stats[:,0].mean())+\ - " Std Dev: "+str(Stats[:,1].mean())+\ - " Good Pts: "+str(Stats[:,2].sum())) - OutTx.write(txVar+"\n") + + print(txVar) + OutTx.write(txVar + "\n") + + OutTx.write("ALL FILES WITH STD DEV < " + str(MINSIG) + " YIELD (" + str(indexGood) + " files)::\n") + txVar = ("Avg Zdr: " + str(Stats[:, 0].mean()) + + " Std Dev: " + str(Stats[:, 1].mean()) + + " Good Pts: " + str(Stats[:, 2].sum())) + OutTx.write(txVar + "\n") OutTx.close() - + # Change the output name of the file - cmd = 'mv zdr_offset.out '+RadName+'_zdr_offset.out' + cmd = 'mv zdr_offset.out ' + RadName + '_zdr_offset.out' os.system(cmd) - + if multiStat: - plot_add_stats(instance='multi', avg=Stats[:,0].mean(),\ - sd=Stats[:,1].mean(), points=Stats[:,2].sum(),\ - good_samples=indexGood, tot_samples=index) + plot_add_stats(instance='multi', avg=Stats[:, 0].mean(), + sd=Stats[:, 1].mean(), points=Stats[:, 2].sum(), + good_samples=indexGood, tot_samples=index) else: - plot_add_stats(instance='single', avg=Stats[0,0],\ - sd=Stats[0,1], points=Stats[0,2]) - + plot_add_stats(instance='single', avg=Stats[0, 0], + sd=Stats[0, 1], points=Stats[0, 2]) + # plt.show() # Create final output file name - fOut = RadName+'_'+pOut+'.'+pType - print "Creating "+ fOut - plt.savefig(fOut,format=pType) + fOut = RadName + '_' + pOut + '.' + pType + print("Creating " + fOut) + plt.savefig(fOut, format=pType) diff --git a/scripts/listfile_raw.py b/scripts/listfile_raw.py index a88e335..78286f9 100755 --- a/scripts/listfile_raw.py +++ b/scripts/listfile_raw.py @@ -14,13 +14,11 @@ 19 Jun 2014 Nick Guy (NRC; NOAA/NSSL) """ -#====================================================================== -# Load the needed packages import pyart.io as fio import argparse import os from glob import glob -#====================================================================== + if __name__ == '__main__': @@ -28,77 +26,74 @@ parser = argparse.ArgumentParser( description='Show file characteristics of raw Sigmet files.') parser.add_argument('searchstring', type=str, help='radar file(s) to list, if more than one file quotations are needed') - parser.add_argument('outFile',type=str, help='output file name to store information') + parser.add_argument('outFile', type=str, help='output file name to store information') igroup = parser.add_argument_group( title='output characteristic file, optional', description=('Decide whether to create output file:')) - + parser.add_argument('-m', '--make_file', action='store_true', help='create an output text file with characteristics') - + args = parser.parse_args() - # ====================================================================== - + # Search for the file(s) flist = glob(args.searchstring) - + # Create an output file if requested if args.make_file: - OutTx = open(args.outFile,'w') - + OutTx = open(args.outFile, 'w') + for file in flist: fLName, fileExt = os.path.splitext(file) fName = os.path.basename(file) - + # Get the file size fSize = os.path.getsize(file) # read in the file - radar = fio.read_sigmet(file,file_field_names=True) + radar = fio.read_sigmet(file, file_field_names=True) - # Create a latitude string + # Create a latitude string if radar.latitude['units'] == 'degrees_north': - rlat = str(radar.latitude['data'])[1:-1]+'N' + rlat = str(radar.latitude['data'])[1:-1] + 'N' if radar.latitude['units'] == 'degrees_south': - rlat = str(radar.latitude['data'])[1:-1]+'S' - + rlat = str(radar.latitude['data'])[1:-1] + 'S' + # Create a longitude string if radar.longitude['units'] == 'degrees_east': - rlon = str(radar.longitude['data'])[1:-1]+'E' + rlon = str(radar.longitude['data'])[1:-1] + 'E' if radar.longitude['units'] == 'degrees_west': - rlon = str(radar.longitude['data'])[1:-1]+'W' - + rlon = str(radar.longitude['data'])[1:-1] + 'W' + # Create field names, check for extended header if radar.metadata['sigmet_extended_header'] == 'true': Vars = "'Xhdr', " + str(radar.fields.keys())[1:-1] else: Vars = str(radar.fields.keys())[1:-1] - + # Create Altitude variable Alt = str(radar.altitude['data'])[1:-1] + ' ' + radar.altitude['units'] - + # Create PRF variable PRF = str(1./radar.instrument_parameters['prt']['data'][0]) + ' Hz' - + # Create Pulse width variable pWid = str(radar.instrument_parameters['pulse_width']['data'][0]) + ' us' - - txVar = (fName+" "+radar.metadata['instrument_name'] + " " + - radar.metadata['sigmet_task_name'] + " File Size: " + str(fSize) + - " Alt: " + Alt + " PRF: " + PRF + - " Pulse Width: " + pWid + " " + Vars + " " + rlat + " " + rlon) - + + txVar = (fName+" "+radar.metadata['instrument_name'] + " " + + radar.metadata['sigmet_task_name'] + " File Size: " + str(fSize) + + " Alt: " + Alt + " PRF: " + PRF + + " Pulse Width: " + pWid + " " + Vars + " " + rlat + " " + rlon) + # Echo this data to the terminal window - print txVar - + print(txVar) + # Add to output file if desired if args.make_file: - OutTx.write(txVar+"\n") - + OutTx.write(txVar + "\n") + del fName, fileExt, radar, rlat, rlon, Vars, Alt, PRF, pWid, txVar - if args.make_file: OutTx.close() - \ No newline at end of file diff --git a/scripts/zdrcal_3panel b/scripts/zdrcal_3panel index 087e49c..886d44a 100755 --- a/scripts/zdrcal_3panel +++ b/scripts/zdrcal_3panel @@ -1,7 +1,5 @@ #! /usr/bin/env python -#====================================================================== -# Import needed libraries import pyart import argparse import numpy as np @@ -10,8 +8,9 @@ import matplotlib.pylab as plt from glob import glob from pyradarmet.zdrcal import calculate_zdr_offset as zdrcal -from pyradarmet.zdrcal import plot_zdr_rhv,plot_add_stats,oplot_zdr_rhv,plot_3panel,oplot_3panel -#====================================================================== +from pyradarmet.zdrcal import (plot_zdr_rhv, plot_add_stats, + oplot_zdr_rhv, plot_3panel, oplot_3panel) + # The user can modify these next variables to control how the program processes data MINSIG = 0.5 # Maximum standard deviation to accept in analysis MINSAMPLE = 100 # Minimum number of samples used to perform calculations @@ -23,8 +22,8 @@ ZH_LIMS = [-20., 30.] # Range of Reflectivity to plot RH_LIMS = [0.8, 1.0] # Range of Copolar Correlation Coefficient to plot DR_LIMS = [-1.5, 1.5] # Range of Differential reflectivity to plot KD_LIMS = [0.2, 0.4] # Range of Specific Differential Phase to plot -PD_LIMS = [-.4,.4] -STD_LIMS = [.1,.7] +PD_LIMS = [-.4, .4] +STD_LIMS = [.1, .7] pOut = "zdr_bias_3panel" pType = "png" @@ -60,27 +59,27 @@ if __name__ == '__main__': parser.add_argument('-v', '--version', action='version', version='Py-ART version %s' % (pyart.__version__)) - + parser.add_argument('-p', '--plot', action='store_true', help='create an output plot file') # parser.add_argument('-t', '--text', action='store_true', # help='create and output text file') - + args = parser.parse_args() # ====================================================================== - + # Search for the file(s) # if args.dir: # os.chdir(args.dir) flist = glob(args.searchstring) - + multiStat = False if len(flist) > 1: multiStat = True # Create an array to hold data, only really used if multiple files - Stats = np.ma.empty([len(flist),3]) + Stats = np.ma.empty([len(flist), 3]) - index = 0 + index = 0 indexGood = 0 for file in flist: fName, fileExt = os.path.splitext(file) @@ -101,25 +100,25 @@ if __name__ == '__main__': radar = pyart.io.read(file) # Call the bias calculation routine - zdr_bias = zdrcal(radar, remove_first_n_gates=GATESKIP, sample_min=MINSAMPLE,\ + zdr_bias = zdrcal(radar, remove_first_n_gates=GATESKIP, sample_min=MINSAMPLE, rhv_min=MINRHV, sig_min=MINSIG, dbz_min=MINDBZ, Htmax=MAXHT) - + # If minimun number of samples is met if zdr_bias is not None: - txVar = ("Avg Zdr: "+str(zdr_bias['volume_average'])+\ - " Std Dev: "+str(zdr_bias['volume_standard_deviation'])+\ - " Good Pts: "+str(zdr_bias['volume_number_good'])+\ - " File: "+os.path.basename(fName)) - + txVar = ("Avg Zdr: " + str(zdr_bias['volume_average']) + + " Std Dev: " + str(zdr_bias['volume_standard_deviation']) + + " Good Pts: " + str(zdr_bias['volume_number_good']) + + " File: " + os.path.basename(fName)) + if zdr_bias['volume_standard_deviation'] > MINSIG: - Stats[index,0] = np.ma.masked - Stats[index,1] = np.ma.masked - Stats[index,2] = np.ma.masked + Stats[index, 0] = np.ma.masked + Stats[index, 1] = np.ma.masked + Stats[index, 2] = np.ma.masked else: - Stats[index,0] = zdr_bias['volume_average'] - Stats[index,1] = zdr_bias['volume_standard_deviation'] - Stats[index,2] = zdr_bias['volume_number_good'] - + Stats[index, 0] = zdr_bias['volume_average'] + Stats[index, 1] = zdr_bias['volume_standard_deviation'] + Stats[index, 2] = zdr_bias['volume_number_good'] + # Create plot if desired if args.plot: if index == 0: @@ -127,43 +126,48 @@ if __name__ == '__main__': #### zdr_bias['height_average'],zdr_bias['kdp_height_average'], ## zdr_bias['height_average'],zdr_bias['height_standard_deviation'], ## zdr_bias['range']['data']/1000.,ymax=MAXHT, - fig, ax1, ax2, ax3 = plot_3panel(zdr_bias['dbz_height_average_filt'], + fig, ax1, ax2, ax3 = plot_3panel( + zdr_bias['dbz_height_average_filt'], ## zdr_bias['height_average_filt'],zdr_bias['phidp_height_average_filt'], - zdr_bias['height_average_filt'],zdr_bias['height_standard_deviation_filt'], - zdr_bias['range']['data']/1000.,ymax=MAXHT, - xlims1=ZH_LIMS,xlims2=DR_LIMS,xlims3=std_lims, - label=True,labelTx=zdr_bias['radar_name'],X3lab='Std Dev') + zdr_bias['height_average_filt'], + zdr_bias['height_standard_deviation_filt'], + zdr_bias['range']['data']/1000., ymax=MAXHT, + xlims1=ZH_LIMS, xlims2=DR_LIMS, xlims3=std_lims, + label=True, labelTx=zdr_bias['radar_name'], X3lab='Std Dev') else: if zdr_bias['height_standard_deviation'].mean() > MINSIG: alff = 0.1 else: alff = 1. - oplot_3panel(zdr_bias['dbz_height_average'],zdr_bias['height_average'], - zdr_bias['height_standard_deviation'],zdr_bias['range']['data']/1000., - ax1=ax1,ax2=ax2,ax3=ax3,alf=alff) + oplot_3panel(zdr_bias['dbz_height_average'], + zdr_bias['height_average'], + zdr_bias['height_standard_deviation'], + zdr_bias['range']['data']/1000., + ax1=ax1, ax2=ax2, ax3=ax3, alf=alff) # Increment indices index += 1 if zdr_bias['height_standard_deviation'].mean() > MINSIG: indexGood += 1 - + RadName = zdr_bias['radar_name'] - + del zdr_bias # If zdr_bias returns None, then too few samples were uncovered. else: txVar = ("Not enough samples in volume") - print txVar + print(txVar) if multiStat: - plot_add_stats(instance='multi',avg=Stats[:,0].mean(), - sd=Stats[:,1].mean(),points=Stats[:,2].sum(), - good_samples=indexGood,tot_samples=index,labelpos=(.4,.83)) + plot_add_stats(instance='multi', avg=Stats[:, 0].mean(), + sd=Stats[:, 1].mean(), points=Stats[:, 2].sum(), + good_samples=indexGood, tot_samples=index, labelpos=(.4, .83)) else: - plot_add_stats(instance='single',avg=Stats[0,0],sd=Stats[0,1],points=Stats[0,2], - labelpos=(.4,.83)) - + plot_add_stats(instance='single', avg=Stats[0, 0], + sd=Stats[0, 1], points=Stats[0, 2], + labelpos=(.4, .83)) + # plt.show() # Create final output file name - fOut = RadName+'_'+pOut+'.'+pType - print "Creating "+ fOut - plt.savefig(fOut,format=pType) + fOut = RadName + '_' + pOut + '.' + pType + print("Creating " + fOut) + plt.savefig(fOut, format=pType) diff --git a/setup.py b/setup.py index 49e83a4..889eb6f 100644 --- a/setup.py +++ b/setup.py @@ -65,58 +65,58 @@ try: basefolder = os.environ["GTOPO_DATA"] - print "Topo data will be installed to %s" % os.environ["GTOPO_DATA"] + print("Topo data will be installed to %s" % os.environ["GTOPO_DATA"]) except KeyError: basefolder = "./pyradarmet/" - print "GTOPO_DATA environmental variable not set, "\ - "using default path in package directory" + print("GTOPO_DATA environmental variable not set, " + "using default path in package directory") archfolder=os.path.join(basefolder, 'gz') datafolder=os.path.join(basefolder, 'data') -names =["antarcps", "e060n40", "e100n40", "e140n40", "w020n40", "w060n90", - "w100n90", "w140n90", "w180s10", "e020n40", "e060n90", "e100n90", - "e140n90", "w020n90", "w060s10", "w100s10", "w140s10", "w180s60", - "e020n90", "e060s10", "e100s10", "e140s10", "w020s10", "w060s60", - "w120s60", "w180n40", "e020s10", "e060s60", "e120s60", "w000s60", - "w060n40", "w100n40", "w140n40", "w180n90"] +names = ["antarcps", "e060n40", "e100n40", "e140n40", "w020n40", "w060n90", + "w100n90", "w140n90", "w180s10", "e020n40", "e060n90", "e100n90", + "e140n90", "w020n90", "w060s10", "w100s10", "w140s10", "w180s60", + "e020n90", "e060s10", "e100s10", "e140s10", "w020s10", "w060s60", + "w120s60", "w180n40", "e020s10", "e060s60", "e120s60", "w000s60", + "w060n40", "w100n40", "w140n40", "w180n90"] for folder in [archfolder, datafolder]: - print folder - if not os.path.exists(folder): - os.makedirs(folder) + print(folder) + if not os.path.exists(folder): + os.makedirs(folder) start_time = time.time() for filename in names: - archfilebase = filename+".tar.gz" - archfile = os.path.join(archfolder, archfilebase) - hdrfilebase = filename+".hdr" - demfilebase = filename+".dem" - - if not os.path.exists(archfile): - print("downloading: {0}".format(archfilebase)) - ftpfile = urllib2.urlopen(gtopourl+archfilebase) - localfile = open(archfile, "wb") - shutil.copyfileobj(ftpfile, localfile) - localfile.close() - - for exfilebase in [hdrfilebase, demfilebase]: - exfile = os.path.join(datafolder,exfilebase) - if not os.path.exists(exfile): - print("extracting: {0} from {1}".format(exfilebase, archfilebase)) - tar = tarfile.open(archfile) - tarobj = tar.extractfile(exfilebase.upper()) - localfile = open(exfile, "wb") - shutil.copyfileobj(tarobj, localfile) - tarobj.close() - localfile.close() - tar.close() + archfilebase = filename+".tar.gz" + archfile = os.path.join(archfolder, archfilebase) + hdrfilebase = filename+".hdr" + demfilebase = filename+".dem" + + if not os.path.exists(archfile): + print("downloading: {0}".format(archfilebase)) + ftpfile = urllib2.urlopen(gtopourl+archfilebase) + localfile = open(archfile, "wb") + shutil.copyfileobj(ftpfile, localfile) + localfile.close() + + for exfilebase in [hdrfilebase, demfilebase]: + exfile = os.path.join(datafolder,exfilebase) + if not os.path.exists(exfile): + print("extracting: {0} from {1}".format(exfilebase, archfilebase)) + tar = tarfile.open(archfile) + tarobj = tar.extractfile(exfilebase.upper()) + localfile = open(exfile, "wb") + shutil.copyfileobj(tarobj, localfile) + tarobj.close() + localfile.close() + tar.close() download_time = time.time() - start_time -print "Total download time: %g seconds"%(download_time) -#print "Removing tar directory: %s" % archfolder -#shutil.rmtree(archfolder) +print("Total download time: %g seconds"%(download_time)) +# print "Removing tar directory: %s" % archfolder +# shutil.rmtree(archfolder) ######### #- Run setup if datafolder == './pyradarmet/data/': @@ -154,7 +154,7 @@ install_requires = ['Numpy >=1.7.2'], ) else: - print "In EXTERNAL" + print("In EXTERNAL") setup( name=PACKNAME, @@ -190,4 +190,4 @@ install_time = time.time() - start_time -print "Total install time: %g seconds"%(install_time) +print("Total install time: %g seconds" % (install_time)) From f8a5d37788445b50697c8dd24d8da1e9a43a0422 Mon Sep 17 00:00:00 2001 From: Nick Guy Date: Wed, 21 Feb 2018 21:57:37 -0800 Subject: [PATCH 2/3] Remove the install script --- .travis.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 2fd4d04..f520663 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,7 +29,6 @@ matrix: before_install: - pip install pytest pytest-cov - pip install coveralls -install: source continuous_integration/install.sh script: eval xvfb-run nosetests $NOSE_ARGS after_success: # Ignore coveralls failures as the coveralls server is not very reliable @@ -37,4 +36,4 @@ after_success: # because the coverage report failed to be published. - if [[ "$COVERALLS" == "true" ]]; then coveralls || echo "failed"; fi # Build docs if requested -- if [[ "$DOC_BUILD" == "true" ]]; then cd $TRAVIS_BUILD_DIR; source continuous_integration/build_docs.sh; fi \ No newline at end of file +- if [[ "$DOC_BUILD" == "true" ]]; then cd $TRAVIS_BUILD_DIR; source continuous_integration/build_docs.sh; fi From 04bb278000ce3565e661aceb53ad282221077217 Mon Sep 17 00:00:00 2001 From: Nick Guy Date: Wed, 21 Feb 2018 22:18:20 -0800 Subject: [PATCH 3/3] Add requirements file --- .travis.yml | 1 + requirements.txt | 11 +++++++++++ 2 files changed, 12 insertions(+) create mode 100644 requirements.txt diff --git a/.travis.yml b/.travis.yml index f520663..f5aae45 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,6 +29,7 @@ matrix: before_install: - pip install pytest pytest-cov - pip install coveralls +- pip install --upgrade -r requirements.txt script: eval xvfb-run nosetests $NOSE_ARGS after_success: # Ignore coveralls failures as the coveralls server is not very reliable diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..9f1d17e --- /dev/null +++ b/requirements.txt @@ -0,0 +1,11 @@ +numpy +matplotlib.cm +glob +os +subprocess +netCDF4 +time +pkg_resources +scipy +mpl_toolkits.basemap +pyart \ No newline at end of file