forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
22 KiB
22 KiB
1.0 Basic Variable Extraction¶
In [ ]:
from __future__ import (absolute_import, division, print_function, unicode_literals)
from wrf import getvar
from netCDF4 import Dataset as nc
ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00")
In [ ]:
p = getvar(ncfile, "P")
print(p)
1.0.1 DataArray attributes: 'dims', 'coords', 'attrs'¶
In [ ]:
print("dims: ", p.dims)
print("coords: ", p.coords)
print("attrs: ", p.attrs)
del p
1.0.2 Removing implicit 'squeeze' behavior to preserve single sized dimensions¶
In [ ]:
p_nosqueeze = getvar(ncfile, "P", timeidx=0, squeeze=False)
print (p_nosqueeze)
1.0.3 Single element metadata¶
In [ ]:
print (p_nosqueeze[0,0,100,200])
del p_nosqueeze
1.0.4 Disabling/Enabling xarray¶
In [ ]:
from wrf import disable_xarray, enable_xarray
# Disable xarray completely
disable_xarray()
p_no_meta = getvar(ncfile, "P")
print(type(p_no_meta))
print (p_no_meta)
del p_no_meta
enable_xarray()
# Disable on extraction
p_no_meta = getvar(ncfile, "P", meta=False)
print("\n")
print(type(p_no_meta))
print(p_no_meta)
del p_no_meta
2.0 Sequences of Input Files¶
2.0.1 Combining via the 'cat' method¶
In [ ]:
wrflist = [ncfile, ncfile, ncfile]
p_cat = getvar(wrflist, "P", method="cat")
print(p_cat)
del p_cat
2.0.2 Combining via the 'join' method¶
In [ ]:
p_join = getvar(wrflist, "P", method="join")
print(p_join)
Note how the Time dimension was replaced with the file dimension, due to the 'squeezing' of the Time dimension.
To maintain the Time dimension, set squeeze to False.
In [ ]:
p_join = getvar(wrflist, "P", timeidx=0, method="join", squeeze=False)
print(p_join)
del p_join
2.0.3 Dictionary Sequences¶
In [ ]:
wrf_dict = {"label1" : [ncfile, ncfile],
"label2" : [ncfile, ncfile]}
p_dict = getvar(wrf_dict, "P")
print(p_dict)
del p_dict
2.0.4 Generator Sequences¶
In [ ]:
def gen_seq():
wrfseq = [ncfile, ncfile, ncfile]
for wrf in wrfseq:
yield wrf
p_gen = getvar(gen_seq(), "P", method="join")
print(p_gen)
del p_gen
2.0.5 Custom Iterable Classes¶
In [ ]:
class FileGen(object):
def __init__(self, ncfile, count=3):
self._total = count
self._i = 0
self.ncfile = [ncfile]*count
def __iter__(self):
return self
def next(self):
if self._i >= self._total:
raise StopIteration
else:
val = self.ncfile[self._i]
self._i += 1
return val
# Python 3
def __next__(self):
return self.next()
obj_gen = FileGen(ncfile, 3)
p_obj_gen = getvar(gen_seq(), "P", method="join", squeeze=False)
print(p_obj_gen)
del p_obj_gen
3.0 WRF Variable Computational Routines¶
In [ ]:
wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz",
"geopt", "helicity", "lat", "lon", "omg", "p", "pressure",
"pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc",
"theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va",
"wa", "uvmet10", "uvmet", "z", "ctt"]
vard = {varname: getvar(ncfile, varname, method="join", squeeze=False) for varname in wrf_vars}
for varname in wrf_vars:
print(vard[varname])
(Note all of the NaNs in the above routines which produce missing values (e.g. cape_2d). xarray always converts all masked_array missing values to NaN in order to work with pandas. To get back the original missing values in a numpy masked_array, you need to use the 'npvalues' method from wrf.)
In [ ]:
from wrf import npvalues
masked_ndarray = npvalues(vard["cape_2d"])
print(masked_ndarray)
del masked_ndarray
In [ ]:
keys = [x for x in vard.keys()]
for key in keys:
del vard[key]
3.1 Interpolation Routines¶
3.1.1 Horizontal Level Interpolation¶
In [ ]:
# 500 MB Heights
from wrf import getvar, interplevel
z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
ht_500mb = interplevel(z, p, 500)
print(ht_500mb)
del ht_500mb, z, p
3.1.2 Vertical Cross Section Interpolation¶
In [ ]:
# Pressure using pivot and angle
from wrf import getvar, vertcross
z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
pivot_point = (z.shape[-2] / 2, z.shape[-1] / 2)
angle = 90.0
p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle)
print(p_vert)
del p_vert
# Pressure using start_point and end_point
start_point = (z.shape[-2]/2, 0)
end_point = (z.shape[-2]/2, -1)
p_vert = vertcross(p, z, start_point=start_point, end_point=end_point)
print(p_vert)
del p_vert, p, z
3.1.3 Interpolate 2D Variable to a Line¶
In [ ]:
# T2 using pivot and angle
from wrf import interpline, getvar
t2 = getvar(ncfile, "T2")
pivot_point = (t2.shape[-2] / 2, t2.shape[-1] / 2)
angle = 90.0
t2_line = interpline(t2, pivot_point=pivot_point, angle=angle)
print(t2_line)
del t2_line
# T2 using start_point and end_point
start_point = (t2.shape[-2]/2, 0)
end_point = (t2.shape[-2]/2, -1)
t2_line = interpline(t2, start_point=start_point, end_point=end_point)
print(t2_line)
del t2_line, t2
3.1.4 Vertical Coordinate Interpolation¶
In [ ]:
from wrf import vinterp, getvar
# Interpolate tk to theta levels
tk = getvar(ncfile, "tk")
interp_levels = [200, 300, 500, 1000]
interp_field = vinterp(ncfile,
field=tk,
vert_coord="theta",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
log_p=True)
print(interp_field)
del interp_field
# Interpolate tk to theta-e levels
interp_levels = [200, 300, 500, 1000]
interp_field = vinterp(ncfile,
field=tk,
vert_coord="eth",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
log_p=True)
print(interp_field)
del interp_field
# Interpolate tk to geopotential height (MSL) levels
interp_levels = [30, 60, 90]
interp_field = vinterp(ncfile,
field=tk,
vert_coord="ght_msl",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
log_p=True)
print(interp_field)
del interp_field
# Interpolate tk to geopotential height (MSL) levels
interp_levels = [30, 60, 90]
interp_field = vinterp(ncfile,
field=tk,
vert_coord="ght_agl",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
log_p=True)
print(interp_field)
del interp_field
# Interpolate tk to pressure levels
interp_levels = [850, 500]
interp_field = vinterp(ncfile,
field=tk,
vert_coord="pressure",
interp_levels=interp_levels,
extrapolate=True,
field_type="tk",
log_p=True)
print(interp_field)
del interp_field, tk
3.2 Lat/Lon to X/Y Routines¶
In [ ]:
from wrf.latlon import get_ll, get_ij # These names are going to change
a = get_ll(ncfile, [400,105], [200,205])
b = get_ij(ncfile, 45.5, -110.8)
# Note: Lists/Dictionaries of files will add a new dimension ('domain') only if the domain is moving
c = get_ll([ncfile, ncfile, ncfile], [400,105], [200,205])
d = get_ll({"label1" : [ncfile, ncfile],
"label2" : [ncfile, ncfile]},
[400,105], [200,205])
print(a)
print("\n")
print(b)
print("\n")
print(c)
print("\n")
print(d)
4.0 Plotting with Cartopy¶
In [ ]:
%matplotlib inline
In [ ]:
# SLP
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from wrf import npvalues, getvar
slp = getvar(ncfile, "slp")
lons = slp.coords["XLONG"]
lats = slp.coords["XLAT"]
wrf_proj = slp.attrs["projection"]
cart_proj = wrf_proj.cartopy()
fig = plt.figure(figsize=(20,20))
ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj)
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
name='admin_1_states_provinces_shp')
ax.add_feature(states, linewidth=.5)
ax.coastlines('50m', linewidth=0.8)
# Can only get this to work if I manually transform the lat/lon points to projected space.
xform_coords = cart_proj.transform_points(crs.PlateCarree(), npvalues(lons), npvalues(lats))
x = xform_coords[:,:,0]
y = xform_coords[:,:,1]
plt.contour(x, y, npvalues(slp), 20, cmap=get_cmap("gist_ncar"))
plt.colorbar(ax=ax, shrink=.7)
ax.set_xlim(wrf_proj.cartopy_xlim())
ax.set_ylim(wrf_proj.cartopy_ylim())
ax.gridlines()
In [ ]:
# 500 MB Heights and Winds
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from wrf import getvar, interplevel, npvalues
p = getvar(ncfile, "pressure")
z = getvar(ncfile, "z", units="dm")
ua = getvar(ncfile, "ua", units="kts")
va = getvar(ncfile, "va", units="kts")
ht_500 = interplevel(z, p, 500)
u_500 = interplevel(ua, p, 500)
v_500 = interplevel(va, p, 500)
lons = ht_500.coords["XLONG"]
lats = ht_500.coords["XLAT"]
wrf_proj = slp.attrs["projection"]
cart_proj = wrf_proj.cartopy()
fig = plt.figure(figsize=(20,20))
ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj)
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
name='admin_1_states_provinces_shp')
ax.add_feature(states, linewidth=0.5)
ax.coastlines('50m', linewidth=0.8)
# Can only get this to work if I manually transform the lat/lon points to projected space.
xform_coords = cart_proj.transform_points(crs.PlateCarree(), npvalues(lons), npvalues(lats))
x = xform_coords[:,:,0]
y = xform_coords[:,:,1]
plt.contour(x, y, npvalues(ht_500), 20, cmap=get_cmap("plasma"))
plt.barbs(x[::50,::50], y[::50,::50], npvalues(u_500[::50, ::50]), npvalues(v_500[::50, ::50]))
plt.colorbar(ax=ax, shrink=.7)
ax.set_xlim(wrf_proj.cartopy_xlim())
ax.set_ylim(wrf_proj.cartopy_ylim())
ax.gridlines()
In [ ]:
# Cross-section of pressure using xarray's builtin plotting
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from wrf import getvar, vertcross, npvalues
p = getvar(ncfile, "pressure")
z = getvar(ncfile, "z", units="dm")
pivot_point = (z.shape[-2] / 2, z.shape[-1] / 2)
angle = 90.0
p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle)
fig = plt.figure(figsize=(20,8))
ax = plt.axes([0.1,0.1,0.8,0.8])
p_vert.plot.contour(ax=ax, levels=[0 + 50*n for n in range(20)], cmap=get_cmap("viridis"))
Multi-time Moving Domain Files¶
In [ ]:
import os
from wrf import getvar
from netCDF4 import Dataset as nc
dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi"
ncfilenames = [os.path.join(dir, x) for x in os.listdir(dir) if x.find("_d02_") > 0]
ncfiles = [nc(x) for x in ncfilenames]
#print (ncfiles[0].variables["XLONG"][0,0,-1], ncfiles[0].variables["XLONG"][-1,0,-1])
#print (ncfiles[1].variables["XLONG"][0,0,-1], ncfiles[1].variables["XLONG"][-1,0,-1])
#print (ncfiles[-1].variables["XLONG"][0,0,-1], ncfiles[-1].variables["XLONG"][-1,0,-1])
In [ ]:
p = getvar(ncfiles, "P", timeidx=3, method="join", meta=True, squeeze=True)
In [ ]:
print (p)
#print (p.attrs["projection"].shape)
In [ ]:
print (p.attrs["projection"])
In [ ]:
ncfiles[2].variables["XTIME"][:]
In [ ]:
p = getvar(ncfiles, "P", timeidx=None, method="cat", meta=True, squeeze=True)
In [ ]:
print (p)
In [ ]:
print (type(p.coords["Time"]))
In [ ]:
import datetime
import pandas
print (type(p.coords["Time"].values.astype(datetime.datetime)))
print (repr(datetime.datetime.utcfromtimestamp(p.coords["Time"][0].values.astype(int) * 1E-9)))
print (pandas.to_datetime(p.coords["Time"].values))
In [ ]:
wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz",
"geopt", "helicity", "lat", "lon", "omg", "p", "pressure",
"pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc",
"theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va",
"wa", "uvmet10", "uvmet", "z", "ctt"]
vard = {varname: getvar(ncfiles, varname, method="join", squeeze=False) for varname in wrf_vars}
for varname in wrf_vars:
print(vard[varname])
In [ ]: