# Copyright (C) 2013-2016 Martin Vejmelka, UC Denver
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
# of the Software, and to permit persons to whom the Software is furnished to do
# so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
# INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR
# A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
import netCDF4
import pytz
from datetime import datetime, timedelta
import numpy as np
import sys
import logging
[docs]class WRFModelData:
"""
This class represents information obtained from a wrfout or wrfinput as generated by WPS or WRF.
Methods for loading data and deriving new variables on the fly are provided.
"""
def __init__(self, path, variables = ['T2', 'Q2', 'PSFC', 'RAINNC', 'RAINC']):
"""
Load data from a wrfinput/wrfout file. The standard loaded fields are those
required for fuel moisture data assimilation, that is: 'T2', 'Q2', 'PSFC', 'RAINNC', 'RAINC'.
In addition, the spatiotemporal information is always loaded: 'XLAT', 'XLONG', 'Times'.
:param path: path to the wrfinput/wrfout
:param variables: variables that should be loaded beyond the standard fields
"""
self.path = path
self.load_data(variables)
[docs] def load_data(self, var_names):
"""
Loads selected variables from the wrfinput/wrfoutput file.
The spatiotemporal information is always loaded: 'XLAT', 'XLONG', 'Times'.
:param var_names: list of variables to load (beyond 'XLAT', 'XLONG' and 'Times')
"""
self.fields = {}
d = netCDF4.Dataset(self.path)
for vname in var_names:
self.fields[vname] = d.variables[vname][:,...]
self.fields['lat'] = d.variables['XLAT'][0,:,:]
self.fields['lon'] = d.variables['XLONG'][0,:,:]
# time is always loaded and encoded as a list of python datetime objects
gmt_tz = pytz.timezone('GMT')
tm = d.variables['Times'][:,...]
tp = []
for t in tm:
dt = datetime.strptime(''.join(t), '%Y-%m-%d_%H:%M:%S')
dt = dt.replace(tzinfo = gmt_tz)
tp.append(dt)
self.fields['GMT'] = tp
d.close()
# if we have all the rain variables, compute the rainfall in each window
if all([v in var_names for v in ['RAINNC', 'RAINC']]):
self.compute_rainfall_per_timestep()
# remove the fields to reduce memory consumption
del self.fields['RAINNC']
del self.fields['RAINC']
else:
# if no rainfall data is available, we set RAIN to zeros
lat = self.fields['lat']
rain_shape = (1, lat.shape[0], lat.shape[1])
self.fields['RAIN'] = np.zeros(rain_shape)
# precompute the equilibrium fields needed everywhere
self.equilibrium_moisture()
[docs] def slice_field(self, field_name):
"""
Remove the temporal dimension from the field by only keeping
the field for the first time instant.
"""
self.fields[field_name] = self.fields[field_name][0,:,:]
[docs] def compute_rainfall_per_timestep(self):
"""
Compute the rainfall per timestep at each grid point from
WRF variables RAINNC and RAINC.
"""
rainnc = self.fields['RAINNC']
rainc = self.fields['RAINC']
rain = np.zeros_like(rainnc)
rain_old = np.zeros_like(rainnc[0,:,:])
tm = self.fields['GMT']
# compute incremental rainfall and store it in mm/hr in the rain variable
for i in range(1, len(tm)):
t1 = tm[i-1]
t2 = tm[i]
dt = (t2 - t1).seconds
rain[i, :, :] = ((rainc[i,:,:] + rainnc[i,:,:]) - rain_old) * 3600.0 / dt
rain_old[:] = rainc[i,:,:]
rain_old += rainnc[i,:,:]
self.fields['RAIN'] = rain
[docs] def get_gmt_times(self):
"""
Returns the local time (depends on time zone set).
"""
return self['GMT']
[docs] def get_lons(self):
"""
Return longitude of grid points.
"""
return self['lon']
[docs] def get_lats(self):
"""
Return lattitute of grid points.
"""
return self['lat']
[docs] def get_field(self, field_name):
"""
Return the field with the name field_name.
"""
return self.field[field_name]
[docs] def get_domain_extent(self):
"""
Return smallest enclosing aligned rectangle of domain.
return is a tuple (min(lon), min(lat), max(lon), max(lat)).
"""
lat = self['lat']
lon = self['lon']
return (np.min(lon), np.min(lat), np.max(lon), np.max(lat))
def __getitem__(self, name):
"""
Access a variable from the fields dictionary.
"""
return self.fields[name]
[docs] def equilibrium_moisture(self):
"""
Uses the fields of the WRF model to compute the equilibrium
field.
"""
# load the standard fields
P = self['PSFC']
Q = self['Q2']
T = self['T2']
self.check_variable(P, 'pressure', 1000, 100000)
self.check_variable(T, 'temperature', 200, 330)
self.check_variable(Q, 'water/vapor ratio', 1e-8, 0.5)
Pi = np.copy(P)
Ti = np.copy(T)
Qi = np.copy(Q)
Pi[1:,:,:] = 0.5 * (P[:-1,:,:] + P[1:,:,:])
Qi[1:,:,:] = 0.5 * (Q[:-1,:,:] + Q[1:,:,:])
Ti[1:,:,:] = 0.5 * (T[:-1,:,:] + T[1:,:,:])
# saturated vapor pressure (at each location, size n x 1)
Pws = np.exp(54.842763 - 6763.22/Ti - 4.210 * np.log(Ti) + 0.000367*Ti + np.tanh(0.0415*(Ti - 218.8))
* (53.878 - 1331.22/Ti - 9.44523 * np.log(Ti) + 0.014025*Ti))
# water vapor pressure (at each location, size n x 1)
Pw = Pi * Qi / (0.622 + (1 - 0.622) * Qi)
# relative humidity (percent, at each location, size n x 1)
H = 100 * Pw / Pws
self.check_variable(H, 'relative humidity', 0., 100.)
mxpos = np.unravel_index(np.argmax(H),H.shape)
logging.info('DIAG: maximum humidity is at %g,%g at time %s.' % (self['lat'][mxpos[1],mxpos[2]],self['lon'][mxpos[1],mxpos[2]], self['GMT'][mxpos[0]]))
H = np.minimum(H, 100.)
# drying/wetting fuel equilibrium moisture contents (location specific,
# n x 1)
d = 0.924*H**0.679 + 0.000499*np.exp(0.1*H) + 0.18*(21.1 + 273.15 - Ti)*(1 - np.exp(-0.115*H))
w = 0.618*H**0.753 + 0.000454*np.exp(0.1*H) + 0.18*(21.1 + 273.15 - Ti)*(1 - np.exp(-0.115*H))
d *= 0.01
w *= 0.01
# this is here to _ensure_ that drying equilibrium is always higher than (or equal to) wetting equilibrium
Ed = np.maximum(d, w)
Ew = np.minimum(d, w)
self.check_variable(Ed, 'drying equilibrium', 0.0, 2.5)
self.check_variable(Ew, 'wetting equilibrium', 0.0, 2.5)
self.fields['Ed'] = Ed
self.fields['Ew'] = Ew
[docs] def get_moisture_equilibria(self):
"""
Return the drying and wetting equilibrium.
"""
return self['Ed'], self['Ew']
[docs] def check_variable(self,V,name,mn,mx):
"""
Check if the variable V is outside the range [mn,mx].
"""
if np.any(V < mn):
pos = np.unravel_index(np.argmin(V), V.shape)
logging.error("Found %s less than %g, min is %g at position %d,%d!" % (name,mn,V[pos],pos[0],pos[1]))
if np.any(V > mx):
pos = np.unravel_index(np.argmax(V), V.shape)
logging.error("Found %s higher than %g, max is %g at position %d,%d!" % (name,mx,V[pos],pos[0],pos[1]))