#!/usr/local/miniconda2/bin/python2.7
#Code: nate_modules.py (Python Module Library)
#Author: Nate Snook (CASA/CAPS/SoM)
#Written: Dec. 2009
#Last modified: Dec. 22, 2009
#---------------------------------------------------------------------------------#
# List of contents:
# grdbas_read(str, netcdf=boolean) -- reads the contents of an ARPS grdbas file (or
# an arps history dump file containing grdbas data and returns grid parameters.
# grdbas_read_patches(str, int, int, netcdf=boolean) -- performs the same function as
# grdbas_read, but for data stored as patches (from mpi runs).
# radarinfo_read(str) -- reads and parses the contents of the ARPS 'radarinfo.dat' file
#---------------------------------------------------------------------------------#
import sys
from math import sqrt
import numpy as np
from pycaps.derive.derive_functions import air_density
from pycaps.util import arps_decompress
from PyNIO import Nio
[docs]def grdbas_read(grdbas_filename, format='hdf'):
"""
Reads in basic grid data from an ARPS grdbas file (in HDF format)
Args:
grdbas_filename: The full path to the grdbas file to read from
format: OPTIONAL -- The format of your input data (currently valid options are 'hdf' (default) and 'netcdf')
Returns:
ctrlat: The latitude of the domain center
ctrlon: The longitude of the domain center
trulat1: The first true latitude value for the lambert conformal map projection
trulat2: The second true latitude value for the lambert conformal map projection
trulon: The true longitude value for the lambert conformal map projection
nx: The number of gridpoints in the east-west direction
ny: The number of gridpoints in the north-south direction
nz: The number of gridpoints in the vertical direction
dx: The grid spacing in the east-west direction
dy: The grid spacing in the north-south direction
width_x: The width of the domain in the east-west direction, in meters
width_y: The width of the domain in the north-south direction, in meters
"""
#Handle special cases first.
#2015 spring experiment CONUS domain, 2-dimensional (y, x) dataset:
if grdbas_filename in ['HWT2015_2D', 'HWT_2015_2D', 'hwt2015_2d', 'hwt_2015_2d']:
ctrlat = 38.50
ctrlon = -97.00
trulat1 = 60.0
trulat2 = 30.0
trulon = -97.0
nx = 1683.
ny = 1155.
nz = 1
dx = 3000.
dy = 3000.
#2015 spring experiment CONUS domain, 3-dimensional (z, y, x) dataset:
elif grdbas_filename in ['HWT2015', 'HWT_2015', 'hwt2015', 'hwt_2015']:
ctrlat = 38.50
ctrlon = -97.00
trulat1 = 60.0
trulat2 = 30.0
trulon = -97.0
nx = 1683.
ny = 1155.
nz = 51
dx = 3000.
dy = 3000.
#If this is not a special case, then open the history file containing grid base data
else:
if format in ['netcdf', 'NetCDF', 'NETCDF']:
grdbasfile = Nio.open_file(grdbas_filename, mode = 'r', options = None, history='', format='netcdf')
grdbas_params = ['CTRLAT','CTRLON','TRUELAT1','TRUELAT2','TRUELON','NX','NY','NZ','DX','DY']
else:
grdbasfile = Nio.open_file(grdbas_filename, mode = 'r', options = None, history='', format='hdf')
grdbas_params = ['ctrlat','ctrlon','trulat1','trulat2','trulon','nx','ny','nz','dx','dy']
#Read in geographic data (center lat, lon; true lat, lon)
ctrlat = float(grdbasfile.attributes[grdbas_params[0]])
ctrlon = float(grdbasfile.attributes[grdbas_params[1]])
trulat1 = float(grdbasfile.attributes[grdbas_params[2]])
trulat2 = float(grdbasfile.attributes[grdbas_params[3]])
trulon = float(grdbasfile.attributes[grdbas_params[4]])
#Read in grid data (nx, ny, nz, dx, dy)
nx = int(grdbasfile.attributes[grdbas_params[5]])
ny = int(grdbasfile.attributes[grdbas_params[6]])
nz = int(grdbasfile.attributes[grdbas_params[7]])
dx = float(grdbasfile.attributes[grdbas_params[8]])
dy = float(grdbasfile.attributes[grdbas_params[9]])
#Calculate derived grid data (width_x, width_y)
width_x = (nx - 1) * dx
width_y = (ny - 1) * dy
Nio.close(grdbasfile)
return ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y
#----------------------------------------------------------------------------------------#
[docs]def grdbas_read_patch(grdbas_filename, xpatches, ypatches, format='hdf'):
"""
Reads in basic grid data from ARPS grdbas patch files (in HDF format)
Args:
grdbas_filename: The full path to the grdbas file to read from
xpatches: The number of patches in the x-direction (in arps.input, nproc_x)
ypatches: The number of patches in the y-direction (in arps.input, nproc_y)
format: OPTIONAL -- The format of your file. Valid options include netcdf and HDF (default HDF).
Returns:
ctrlat: The latitude of the domain center
ctrlon: The longitude of the domain center
trulat1: The first true latitude value for the lambert conformal map projection
trulat2: The second true latitude value for the lambert conformal map projection
trulon: The true longitude value for the lambert conformal map projection
nx: The number of gridpoints in the east-west direction in the full domain
ny: The number of gridpoints in the north-south direction in the full domain
nz: The number of gridpoints in the vertical direction in the full domain
dx: The grid spacing in the east-west direction
dy: The grid spacing in the north-south direction
width_x: The width of the domain in the east-west direction, in meters
width_y: The width of the domain in the north-south direction, in meters
nx_patch: The number of gridpoints in the east-west direction in a single patch
ny_patch: The number of gridpoints in the north-south direction in a single patch
"""
#Open the history file containing grid base data
grdbas_patchname = grdbas_filename + '_001001'
if arguments.format in ['netcdf', 'NetCDF', 'NETCDF']:
grdbasfile = Nio.open_file(grdbas_patchname, mode = 'r', options = None, history='', format='netcdf')
grdbas_params = ['CTRLAT', 'CTRLON', 'TRUELAT1', 'TRUELAT2', 'TRUELON', 'NX', 'NY', 'NZ', 'DX', 'DY']
else:
grdbasfile = Nio.open_file(grdbas_patchname, mode = 'r', options = None, history='', format='hdf')
grdbas_params = ['ctrlat', 'ctrlon', 'trulat1', 'trulat2', 'trulon', 'nx', 'ny', 'nz', 'dx', 'dy']
#Read in geographic data (center lat, lon; true lat, lon)
ctrlat = float(grdbasfile.attributes[grdbas_params[0]])
ctrlon = float(grdbasfile.attributes[grdbas_params[1]])
trulat1 = float(grdbasfile.attributes[grdbas_params[2]])
trulat2 = float(grdbasfile.attributes[grdbas_params[3]])
trulon = float(grdbasfile.attributes[grdbas_params[4]])
#Read in grid data (nx, ny, nz, dx, dy)
nx_patch = int(grdbasfile.attributes[grdbas_params[5]])
ny_patch = int(grdbasfile.attributes[grdbas_params[6]])
nz = int(grdbasfile.attributes[grdbas_params[7]])
dx = float(grdbasfile.attributes[grdbas_params[8]])
dy = float(grdbasfile.attributes[grdbas_params[9]])
nx = ((nx_patch - 3) * xpatches) + 3
ny = ((ny_patch - 3) * ypatches) + 3
#Calculate derived grid data (width_x, width_y)
width_x = (nx - 1) * dx
width_y = (ny - 1) * dy
Nio.close(grdbasfile)
return ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y, nx_patch, ny_patch
#----------------------------------------------------------------#
[docs]def read_xy_slice(field, source, level, format='hdf', **kwargs):
"""
Reads the data needed for xyplot to generate an x-y variable field plot.
Args:
field: The variable name associated with the field to be plotted (e.g. u, v, pt, qr)
source: The FULL PATH to the file containing the data to be plotted.
level: The vertical (k) model layer for which data should be read and plotted.
format: OPTIONAL -- The input data format (currently valid options are 'hdf' (default) and 'netcdf')
grdbas: OPTIONAL -- A history file to read grid information from (if different from source)
h2: OPTIONAL -- 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).
truref: OPTIONAL -- 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.
decompress: OPTIONAL -- 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:
var_sfc: The 2D slice of data to be plotted
"""
#----------------------------------------------------#
# EXTRA SHAPEFILE ID LIST #
#----------------------------------------------------#
extra_shape = []
shapevar = []
shapefile_dir = '/data5/nsnook/python_scripts/shapefiles/'
verbose = kwargs.get('verbose', True)
#Determine the appropriate format for reading in data.
if 'format' in kwargs:
if (kwargs['format'] == 'netcdf' or kwargs['format'] =='NetCDF'):
filefmt = 'netcdf'
elif (kwargs['format'] == 'hdf' or kwargs['format'] =='HDF'):
filefmt = 'hdf'
else:
sys.exit('Invalid format error! The only valid formats are hdf and netcdf.')
else:
filefmt = 'hdf'
#Read grid data using the grdbas_read function:
if 'grdbas' in kwargs:
ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y = grdbas_read(kwargs['grdbas'], format=filefmt)
else:
try:
ctrlat, ctrlon, trulat1, trulat2, trulon, nx, ny, nz, dx, dy, width_x, width_y = grdbas_read(source, format=filefmt)
except:
sys.exit('FATAL ERROR: No grdbas file given, and grid data not available in history data!')
#----------------------------------------------------#
# OPEN ARPS HISTORY FILE AND READ DATA #
#----------------------------------------------------#
#Initialize array to store data for plotting
var_sfc_filtered = np.zeros((nx,ny))
dumpfile = Nio.open_file(source, mode = 'r', options = None, history='', format=filefmt)
if 'h2' in kwargs:
dumpfile2 = Nio.open_file(kwargs['h2'], mode = 'r', options = None, history='', format=filefmt)
if (field == 'tsoil' or field == 'qsoil'):
if verbose:
print 'Special case -- Soil variables (tsoil and qsoil) detected.'
var_4D = dumpfile.variables[field][:,:,:,:]
var_3D = var_4D[0,:,:,:]
#If the variable is directly in the file, just get it from there.
if field in dumpfile.variables.keys():
#Read in the data for qr(k,j,i)
var_sfc = dumpfile.variables[field][level,:,:]
if 'decompress' in kwargs:
if verbose:
print "Decompressing ARPS compressed data..."
var_sfc = arps_decompress(dumpfile.variables[field])
var_sfc = var_sfc[level,:,:]
else:
if verbose:
print 'Variable ' + str(field) + ' not found in data. Trying derived variables.'
if (field == 'ptprt') or (field == 'ptpert'): # If the user wants perturbation potential temperat ure (ptprt)...
pt = dumpfile.variables['pt'][level,:,:] #...then calculate it from potential temperature ( pt) and ptbar.
ptbar = dumpfile.variables['ptbar'][level,:,:]
if 'decompress' in kwargs:
pt = arps_decompress(dumpfile.variables['pt'])
pt = pt[level,:,:]
ptbar = arps_decompress(dumpfile.variables['ptbar'])
ptbar = ptbar[level,:,:]
ptprt = np.zeros(ptbar.shape)
ptprt = pt - ptbar #perturbation = actual value - average value
var_sfc = ptprt
elif (field == 'AGL'):
sfc_height = dumpfile.variables['zpsoil'][0,:,:] #...calculate it from zpsoil[sfc] and zp.
zp = dumpfile.variables['zp'][level,:,:]
if 'decompress' in kwargs:
sfc_height = arps_decompress(dumpfile.variables['zpsoil'])
sfc_height = sfc_height[0,:,:]
zp = arps_decompress(dumpfile.variables['zp'])
zp = zp[level,:,:]
agl = zp - sfc_height
var_sfc = agl
elif field in ['mh', 'mg', 'mr', 'mi', 'ms', 'mc']:
p = dumpfile.variables['p'][level,:,:]
pt = dumpfile.variables['pt'][level,:,:]
if 'decompress' in kwargs:
p = arps_decompress(dumpfile.variables['p'])
p = p[level,:,:]
pt = arps_decompress(dumpfile.variables['pt'])
pt = pt[level,:,:]
rho_air = air_density(p = p, pt = pt)
#Determine appropriate variable to read:
if field == 'mh':
qstr = 'qh'
elif field == 'mg':
qstr = 'qg'
elif field == 'mr':
qstr = 'qr'
elif field == 'mi':
qstr = 'qi'
elif field == 'ms':
qstr = 'qs'
elif field == 'mc':
qstr = 'qc'
else:
print str(field) + ' is not a valid microphysical species!'
sys.exit('ERROR: INVALID MICROPHYSICAL SPECIES CHOSEN')
#Read mixing ratio and compute hydrometeor mass:
qx = dumpfile.variables[qstr][level,:,:]
if 'decompress' in kwargs:
qx = arps_decompress(dumpfile.variables[qstr])
qx = qx[level,:,:]
mx = qx * rho_air
var_sfc = 1000 * mx
elif (field == 'Z'): # Lin reflectivity
if 'truref' in kwargs:
if verbose:
print 'using truref from ossedata...'
truref_file = Nio.open_file(kwargs['truref'], mode = 'r', options = None, history='', format='hdf')
var_sfc = truref_file.variables['truref'][level,:,:]
elif 'decompress' in kwargs:
pt = arps_decompress(dumpfile.variables['pt'])
pt = pt[level,:,:] # Potential temperature (K)
p = arps_decompress(dumpfile.variables['p'])
p = p[level,:,:] # Pressure
qr = arps_decompress(dumpfile.variables['qr'])
qr = qr[level,:,:] # Rain water mixing ratio
qs = arps_decompress(dumpfile.variables['qs'])
qs = qs[level,:,:] # Snow mixing ratio
qh = arps_decompress(dumpfile.variables['qh'])
qh = qh[level,:,:] # Hail mixing ratio
n0rain = float(dumpfile.attributes['n0rain'])
n0snow = float(dumpfile.attributes['n0snow'])
n0hail = float(dumpfile.attributes['n0hail'])
rhoice = float(dumpfile.attributes['rhoice'])
rhosnow = float(dumpfile.attributes['rhosnow'])
rhogrpl = float(dumpfile.attributes['rhogrpl'])
rhohail = float(dumpfile.attributes['rhohail'])
linrefl = reflectivity_lin(pt=pt,p=p,qr=qr,qs=qs,qh=qh,n0rain=n0rain,n0snow=n0snow,n0hail=n0hail,rhoice=rhoice,rhosnow=rhosnow,rhogrpl=rhogrpl,rhohail=rhohail)
var_sfc = linrefl
else:
pt = dumpfile.variables['pt'][level,:,:] # Potential temperature (K)
p = dumpfile.variables['p'][level,:,:] # Pressure
qr = dumpfile.variables['qr'][level,:,:] # Rain water mixing ratio
qs = dumpfile.variables['qs'][level,:,:] # Snow mixing ratio
qh = dumpfile.variables['qh'][level,:,:] # Hail mixing ratio
n0rain = float(dumpfile.attributes['n0rain'])
n0snow = float(dumpfile.attributes['n0snow'])
n0hail = float(dumpfile.attributes['n0hail'])
rhoice = float(dumpfile.attributes['rhoice'])
rhosnow = float(dumpfile.attributes['rhosnow'])
rhogrpl = float(dumpfile.attributes['rhogrpl'])
rhohail = float(dumpfile.attributes['rhohail'])
linrefl = reflectivity_lin(pt=pt,p=p,qr=qr,qs=qs,qh=qh,n0rain=n0rain,n0snow=n0snow,n0hail=n0hail,rhoice=rhoice,rhosnow=rhosnow,rhogrpl=rhogrpl,rhohail=rhohail)
linrefl = reflectivity_lin(pt=pt,p=p,qr=qr,qs=qs,qh=qh,n0rain=n0rain,n0snow=n0snow,n0hail=n0hail,rhoice=rhoice,rhosnow=rhosnow,rhogrpl=rhogrpl,rhohail=rhohail)
linrefl = reflectivity_lin(pt=pt,p=p,qr=qr,qs=qs,qh=qh,n0rain=n0rain,n0snow=n0snow,n0hail=n0hail,rhoice=rhoice,rhosnow=rhosnow,rhogrpl=rhogrpl,rhohail=rhohail)
var_sfc = linrefl
elif (field == 'Zdm' or field == 'Z2' or field == 'Z2m'):
if 'truref' in kwargs:
if verbose:
print 'using truref from ossedata...'
truref_file = Nio.open_file(kwargs['truref'], mode = 'r', options = None, history='', format='hdf')
var_sfc = truref_file.variables['truref'][level,:,:]
else:
if verbose:
print 'Variable is dual-moment reflectivity (MY2). Calculating using dropsizedist.py...'
p = dumpfile.variables['p'][level]
pt = dumpfile.variables['pt'][level]
temperature = theta2Temperature(p=p, pt=pt)
density_air = p / (temperature * 287.)
psdc_dm = PSDCollection(
rain=RainPSD(dumpfile.variables['qr'][level], num_concentration=dumpfile.variables['nr'][level]),
snow=SnowPSD(dumpfile.variables['qs'][level], num_concentration=dumpfile.variables['ns'][level]),
hail=HailPSD(dumpfile.variables['qh'][level], num_concentration=dumpfile.variables['nh'][level]),
graupel=GraupelPSD(dumpfile.variables['qg'][level], num_concentration=dumpfile.variables['ng'][level]))
var_sfc = psdc_dm.reflectivity(density_air, method='jzx')
elif (field == 'divergence' or field == 'div'):
if verbose:
print 'Calculating horizontal divergenge...'
u = dumpfile.variables['u'][level,:,:] # u-component of velocity
v = dumpfile.variables['v'][level,:,:] # v-component of velocity
div_sfc = np.zeros(u.shape)
div_sfc[1:-1,1:-1] = ((u[2:,1:-1] - u[:-2,1:-1])/(2 * dx)) + ((v[1:-1,2:] - v[1:-1,:-2])/(2 * dy))
var_sfc = div_sfc
elif (field == 'wspd' or field == 'windspeed'):
u = dumpfile.variables['u'][level,:,:]
v = dumpfile.variables['v'][level,:,:]
var_sfc = sqrt((u * u) + (v * v))
elif arguments.var in ['uh', 'UH', 'updraft_helicity']:
if verbose:
print 'Calculating updraft helicity...'
u = dumpfile.variables['u'][:,:,:]
v = dumpfile.variables['v'][:,:,:]
w = dumpfile.variables['w'][:,:,:]
zp = dumpfile.variables['zp'][:,:,:]
var_sfc = updraft_helicity(u, v, w, dx, dy, zp)
else:
try:
var_sfc = dumpfile.variables[arguments.var][:,:]
except:
print 'Variable ' + str(arguments.var) + ' does not conform to expected dimensions (2D or 3D)'
print ' '
sys.exit('ERROR: Invalid variable shape!')
return var_sfc
#------------------------------------------------------------------------#