#!/usr/bin/env python2.7
#Code: plot_arps_xy.py (Python Script)
# (taken from plot_var_advanced.py v2.1)
#Author: Nate Snook (CASA/CAPS/SoM)
#Date Created: Mar. 2008
#Modification history:
# 20 April 2012 by Robin Tanamachi (Various edits)
# 26 Feb. 2013 by Nate Snook (Improved command-line input and parsing)
# 26 Nov. 2013 by Nate Snook (Updated script to use argparse for command-line input)
# 19 Feb. 2014 by Nate Snook (Added support for NetCDF files)
# 11 Apr. 2014 by Nate Snook (Added support for overlaying extra shapefiles and/or vertical vorticity)
# 11 Apr. 2014 by Nate Snook (Added subdomain plotting options on the command line)
# 13 May 2014 by Nate Snook (Added command-line option to choose the colormap)
# 14 May 2014 by Nate Snook (Fixed bugs in colormap implementation -- thanks to Jon Labriola for bug report)
# 7 Oct. 2014 by Nate Snook (3D-->2D read-in for efficiency, 'from numpy import *' --> 'import numpy as np')
# 21 Apr. 2015 -- Version 2.0 update by Nate Snook:
# Moved colormap specification to a module in ns_modules.py
# Added pcolormesh support through '--pixels' option
# The '--grdbas' argument is now required only if grid data missing in arps file.
# 9 Jun. 2015 by Nate Snook (Added support for compressed ARPS variables (compression opt. 5))
# 14 Jul. 2015 by Nate Snook (Added code to plot an arbitrary line ("--line_at_y" option))
# 1 Oct. 2015 -- Version 2.1 update by Nate Snook:
# Code now uses 'grdbas_read' from 'ns_modules.py' instead of local code to read grdbas data.
# Configuration script 'plot_config.py' is NO LONGER SUPPORTED; must use command-line.
# Updates to comments and in-code documentation.
# 13 Oct. 2015 by Nate Snook (Fixed broken shapefile directory caused by v2.1 update)
# 20 Oct. 2015 by Nate Snook (converted code to plot_arps_xy.py for use in pycaps
# 3 Nov. 2015 by Nate Snook (Converted stand-alone script into a function for use in pycaps)
#
[docs]def xyplot(data, field, source, level, format='hdf', **kwargs):
"""
Plots a single HDF or NetCDF variable field (or field derived from history data)
in the x-y plane, with options for various overlays. Grid data, such as that stored in
an ARPS grdbas file, are required. Lambert conformal map projection will be used.
Args:
data: The 2D x-y slice to be plotted (e.g. u, v, pt, qr). Can be generated quickly using read_xy_slice.
field: The variable name (e.g. 'u', 'v', 'pt', 'qr'). Needed to properly format your plot. Enclose it in single quotes.
source: The FULL PATH to the file containing the data for plotting.
level: The vertical (k) model layer on which to plot an x-y slice.
format: OPTIONAL -- The file format for your input (currently, valid choices are "hdf" (default) and "netcdf")
**kwargs:
grdbas: The history file from which to read grid information (if different from 'source')
h2: A second source file containing other data needed for plotting, if you have one. (example: Read reflectivity from one file, and wind data from another).
name: A name to identify the data (may be used in plot title and output filename)
truref: If reading from an ARPS file processed using ossedata (affects variable name and filename) this should be set to the full path to the truref file.
addshape: Specify a shapefile to overlay data from. Requires the shapefile (and options) to be included in the shapefiles block in plot_arps_xy.py within the pycaps package.
colormap: Used to specify a colormap from ns_colormap.py for use with this plot. If not provided, a default colormap choice appropriate to the variable will be used.
pixels: A flag to set to True if you want a pixel (colormesh) plot instead of a contour plot.
tornado: A flag to set to True if you are plotting values from a high-resolution tornadic dataset. Will choose appropriate values for relevant variables (e.g. wspd, p)
imin: The i-coordinate of the west edge of your plot. If not given, default = 0
imax: The i-coordinate of the east edge of your plot. If not given, default = nx
jmin: The j-coordinate of the south edge of your plot. If not given, default = 0
jmax: The j-coordinate of the north edge of your plot. If not given, default = ny
qcbelow: If given, values below this threshold will be marked as missing/bad data.
qcabove: If given, values above this threshold will be marked as missing/bad data.
urban: A flag to set to True if you want to plot urban data (if not given, False by default).
counties: A flag to set to False if you don't want to plot counties (if not given, True by default).
states: A flag to set to False if you don't want to plot states (if not given, True by default).
coastlines: A flag to set to False if you don't want to plot coastlines (if not given, True by default).
lat: If given, will plot lines of latitude at the given interval (in degrees).
lon: If given, will plot lines of longitude at the given interval (in degrees).
decompress: A flag to set to True if you're using ARPS data with 32 bit integers mapped to 16 bit integers. Will call the decompression function if set to True.
Returns:
<<nothing>> (Generates and saves the requested plot to a .png file)
"""
import matplotlib.pyplot as plt
import sys
import os
from math import sqrt
import numpy as np
from mpl_toolkits.basemap import Basemap
from PyNIO import Nio
#Replaced custom module calls with their pycaps equivalents.
#from ns_modules import air_density, from_levels_and_colors, thresh_setup, grdbas_read
from pycaps.derive import air_density
from pycaps.io import grdbas_read
from pycaps.plot import from_levels_and_colors, thresh_setup #or use matplotlib's from_levels_and_colors
from pycaps.plot.ns_colormap import *
from pycaps.util import arps_decompress
#Read grid data using the grdbas_read function:
verbose = kwargs.get('verbose', True)
if 'grdbas' in kwargs:
if verbose:
print 'attempting to read grdbas file ' + str(kwargs['grdbas']) + '...'
ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y = grdbas_read(kwargs['grdbas'], format=format)
else:
try:
if verbose:
print 'attempting to read grdbas data from source file ' + str(source) + '...'
ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y = grdbas_read(source, format=format)
except:
try:
ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y = grdbas_read(source)
except:
sys.exit('FATAL ERROR: grdbas file not given, and grid data not available in history data!')
#Determine experiment name to use for plotting
if 'name' in kwargs:
if verbose:
print 'using specified experiment name: ' + kwargs['name']
title_filename = kwargs['name']
output_filename = str(kwargs['name']) + '_' + source[-6:] + '_' + str(field) + '_lev%03d.png'%int(level)
else:
if verbose:
print 'No experiment name specified -- using full experiment path instead.'
title_filename = source.split('/')[-2] + '_' + source.split('/')[-1]
output_filename = source.split('/')[-2] + '_' + source.split('/')[-1]+ '_' + str(field) + '_lev%03d.png'%int(level)
if verbose:
print '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-='
print 'Plotting at grid level ' + str(int(level) + 1) + '.'
print '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-='
#Filter bad data...
if 'qcabove' in kwargs:
data = np.where(data > kwargs['qcabove'], -99.0, data)
if 'qcbelow' in kwargs:
data = np.where(data < kwargs['qcbelow'], -99.0, data)
#Take care of missing data for gridtilt reflectivity, Zdr, rho-hv:
if (field in ["outgridtilt_zhh", "outgridtilt_zdr", "outgridtilt_rhv", "obsordr_zhh", "obsordr_zdr", "obsordr_rhv"]):
if verbose:
print 'Ensuring proper treatment of missing values in gridtilt radar field...'
data = np.where(data > -99.0, data, 0.0)
if verbose:
print 'Plotting ' + str(field)
# Replace with nanmax and nanmin, so that NaNs don't bollux up the threshold calculations.
print 'Maximum value = ' + str(np.nanmax(data)) + ' || Minimum value = ' + str(np.nanmin(data))
#Create a figure
x = np.arange(0,width_x+dx,dx) #the '+dx' or '+dy' is so that np.arange() will include the final...
y = np.arange(0,width_y+dy,dy) #...gridpoint at width_x and width_y. Trust me, it's necessary.
x, y = np.meshgrid(x,y)
#lambert conformal map projection containing states, counties, rivers, and urban boundaries
map = Basemap(projection='lcc', width=width_x, height=width_y, lat_1=trulat1, lat_2=trulat2, lat_0=ctrlat, lon_0=ctrlon, resolution='h', area_thresh=1000.) #lambert conformal -- uses parameters from grdbas file!
#----------------------------------------------------#
# SUBDOMAIN SPECIFICATION #
#----------------------------------------------------#
if (('imin' in kwargs) or ('imax' in kwargs) or ('jmin' in kwargs) or ('jmax' in kwargs)):
if 'imin' in kwargs:
xmin = kwargs['imin']
else:
xmin = 0
if 'jmin' in kwargs:
ymin = kwargs['jmin']
else:
ymin = 0
if 'imax' in kwargs:
xmax = kwargs['imax']
else:
xmax = nx
if 'jmax' in kwargs:
ymax = kwargs['jmax']
else:
ymax = ny
#Determine center lat, lon of the subdomain:
xctr = int(((xmax + xmin) / 2)*dx)
yctr = int(((ymax + ymin) / 2)*dy)
lonctr, latctr = map(xctr, yctr, inverse = True)
#Trim arrays (including overlays, if applicable)
data = data[ymin:ymax, xmin:xmax]
if (('qc_below' in kwargs) or ('qc_above' in kwargs)):
data_filtered = data_filtered[ymin:ymax, xmin:xmax]
if 'vortoverlay' in kwargs:
vort_sfc = vort_sfc[ymin:ymax, xmin:xmax]
if 'wind' in kwargs:
u_sfc_kts = u_sfc_kts[ymin:ymax, xmin:xmax]
v_sfc_kts = v_sfc_kts[ymin:ymax, xmin:xmax]
width_x = float(xmax - xmin - 1)*dx
width_y = float(ymax - ymin - 1)*dy
x = np.arange(0, width_x+dx, dx)
y = np.arange(0, width_y+dy, dy)
x, y = np.meshgrid(x,y)
if verbose:
print "----- ----- -----"
print 'Plot will show the following subdomain:'
print 'x: (' + str(xmin) + ',' + str(xmax) + ')'
print 'y: (' + str(ymin) + ',' + str(ymax) + ')'
print 'New ctrlat, ctrlon: (' + str(latctr) + ',' + str(lonctr) + ')'
print "----- ----- -----"
map = Basemap(projection='lcc', width=width_x, height=width_y, lat_1 = trulat1, lat_2 = trulat2, lat_0 = latctr, lon_0 = lonctr, resolution='h', area_thresh=10.) #New lambert conformal map for subdomain.
#Now, handle the actual figure plotting.
fig=plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
if 'tornado' in kwargs:
tornado = kwargs['tornado']
else:
tornado = False
#Set up thresholds, colormap, and colorbar (using thresh_setup in ns_modules.py):
thresholds, cb_ticks, colormap, no_proportional = thresh_setup(field, data.max(), data.min(), tornado)
#Handle masking for those variables that require it:
if field in ['zdr', 'zdr_3d', 'zdr_2d', 'Zdr']: # Differential reflectivity
#For Zdr, mask missing values (typically -999, but could be -99)
data = np.ma.masked_where(data < -98., data)
if field in ['rhohv', 'rho_hv', 'rhv_2d', 'rhv_3d']:
#Cross-polarization correlation coefficient (rho-hv) -- missing value = -1.0
data = np.ma.masked_where(data < -0.99, data)
if field in ['vr', 'radv2d', 'radv3d', 'radvel']:
#For radial velocity, mask missing values (typically -999, but could be -99)
data = np.ma.masked_where(data < -98., data)
#Special setup for mean-mass diameter variables (convert m to mm)
if field in ['mmdi_h', 'mmdi_g', 'mmdi_r', 'mmdi_s', 'mmdi_i', 'mmdi_c']:
data = data * 1000. #m --> mm
#If the user specified a colormap on the command-line, use that one.
if 'colormap' in kwargs:
colormap = colormap_list[kwargs['colormap']]
#create a contour plot
if 'pixels' in kwargs:
if verbose:
print 'Plotting pcolormesh (pixel) map:'
print ' colormap = ' + str(colormap)
print ' thresholds = ' + str(thresholds)
cmap, norm = from_levels_and_colors(thresholds, colormap, extend='neither')
cmap.set_under(colormap[0])
cmap.set_over(colormap[-1])
CF = map.pcolormesh(x, y, data, cmap=cmap, norm=norm)
else:
if verbose:
print 'Now drawing contour plot...'
print ' Colors = ' + str(colormap)
CF = map.contourf(x, y, data, thresholds, colors=colormap, extend='both')
if no_proportional == True:
CB = plt.colorbar(CF, shrink=0.8, extend='both', ticks=cb_ticks)
else:
CB = plt.colorbar(CF, shrink=0.8, extend='both', spacing='proportional', ticks=cb_ticks)
if (field == "Z"):
if verbose:
print '...entering colormap specifications for Z...'
CB.set_label("Radar Reflectivity (dBZ)")
CF.cmap.set_under('#FFFFFF')
CF.cmap.set_over('#BB55DD')
elif field in ['rho-hv', 'rho_hv', 'rhv_2d', 'rhv_3d']:
CB.set_label("Cross-polarization Correlation Coefficient")
CF.cmap.set_under('#FFFFFF')
elif field in ['Zdr', 'ZDR', 'zdr_2d', 'zdr_3d']:
CB.set_label("Differential Reflectivity (dB)")
elif (field == "pt"):
CB.set_label("Potential Temperature (K)")
elif (field == "zp"):
CB.set_label("Meters above Sea Level")
elif (field == "u" or field == "v" or field == "w"):
CB.set_label("m s$^{-1}$")
elif (field == 'vort' or field == "xvort" or field == "yvort" or field == "zvort" or field == "strvort"):
CB.set_label("x 10$^{-5}$ s$^{-1}$")
elif (field == "solx" or field == "soly" or field == "solstr" or field == "tilt" or field == "stretch"):
CB.set_label("x 10$^{-5}$ s$^{-2}$")
elif (field == "sol" or field == "solz"):
CB.set_label("x 10$^{-7}$ s$^{-2}$")
shapefile_dir = os.path.join(os.path.dirname(__file__), '..', 'data')
if 'addshape' in kwargs:
if verbose:
print 'Plotting extra shapefile overlay(s)...'
if (kwargs['addshape'] == '20may2013' or kwargs['addshape'] == '130520'):
map.readshapefile(shapefile_dir + '/20may2013/extractDamagePolys', 'tornado_track', drawbounds= False)
xx, yy = zip(*map.tornado_track[2])
map.plot(xx,yy,linewidth=1.3, color='#000000')
else:
for index, shapefile in enumerate(extra_shape):
map.readshapefile(extra_shape[index], shapevar[index], drawbounds=True, linewidth=0.4, color ='#553333')
if 'counties' in kwargs:
if kwargs['counties'] == False:
pass
else:
if verbose:
print 'reading counties from shapefile'
map.readshapefile(shapefile_dir + '/counties/countyp020','counties',drawbounds=True, linewidth=0.8, color='#333333') #Draws US county boundaries.
else:
if verbose:
print 'reading counties from shapefile'
map.readshapefile(shapefile_dir + '/counties/countyp020','counties',drawbounds=True, linewidth=0.8, color='#333333') #Draws US county boundaries.
if 'states' in kwargs:
if kwargs['states'] == 'False':
pass
else:
map.drawstates(linewidth=1.0, color='#222222')
else:
map.drawstates(linewidth=1.0, color='#222222')
if 'urban' in kwargs:
if kwargs['urban'] == True:
if verbose:
print 'reading urban boundaries from shapefile'
map.readshapefile(shapefile_dir + '/urban/cb_2012_us_uac10_500k', 'urban_areas', drawbounds=True, linewidth=0.5, color='purple') #Draws urban areas in purple.
if 'hwys' in kwargs:
if kwargs['hwys'] == True:
if verbose:
print 'reading major highways from shapefile'
map.readshapefile(shapefile_dir + '/highway/road_l', 'major_highways', drawbounds=True, linewidth=0.7, color='#0000AA')
if 'countries' in kwargs:
if kwargs['countries'] == False:
pass
else:
map.drawcountries(linewidth=1.5, color='#000000')
else:
map.drawcountries(linewidth=1.5, color='#000000')
if 'coastlines' in kwargs:
if kwargs['coastlines'] == False:
pass
else:
map.drawcoastlines(linewidth=2.0, color='#000000')
else:
map.drawcoastlines(linewidth=2.0, color='#000000')
#title the plot
plt.title('Plot of ' + str(field) + ' in ' + str(title_filename) + ' -- Grid level : ' + str(level))
map.drawmapboundary()
#If desired, draw lines for longitude and latitude
if 'lon' in kwargs:
meridians = np.arange(-160., 160., kwargs['lon']) #List of meridians for drawing.
map.drawmeridians(meridians, labels=[0,0,0,1], linewidth=1.0) #Draws meridians (lines of constant longitude)
if 'lat' in kwargs:
parallels = np.arange(0., 80., kwargs['lat']) #List of parallels for drawing.
map.drawparallels(parallels, labels=[1,0,0,0], linewidth=1.0) #Draws parallels (lines of constant latitude)
if verbose:
#Show figure or save as .png:
print "Plot saved to " + str(output_filename)
plt.savefig(output_filename, figsize = (13, 13), dpi=300)
#end if
# Close open HDF files
try:
dumpfile.close()
except:
pass