<jupyter_start><jupyter_text>1.0 Basic Variable Extraction<jupyter_code>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")
ncfile = nc("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00")

p = getvar(ncfile, "P")
print(p)<jupyter_output><empty_output><jupyter_text>1.0.1 DataArray attributes:  'dims', 'coords', 'attrs'<jupyter_code>print("dims: ", p.dims)
print("coords: ", p.coords) 
print("attrs: ", p.attrs)
del p<jupyter_output><empty_output><jupyter_text>1.0.2 Removing implicit 'squeeze' behavior to preserve single sized dimensions<jupyter_code>p_nosqueeze = getvar(ncfile, "P", timeidx=0, squeeze=False)
print (p_nosqueeze)<jupyter_output><empty_output><jupyter_text>1.0.3 Single element metadata<jupyter_code>print (p_nosqueeze[0,0,100,200])
del p_nosqueeze<jupyter_output><empty_output><jupyter_text>1.0.4 Disabling/Enabling xarray<jupyter_code>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<jupyter_output><empty_output><jupyter_text>2.0 Sequences of Input Files   2.0.1 Combining via the 'cat' method<jupyter_code>from wrf import ALL_TIMES
import numpy as np
wrflist = [ncfile, ncfile, ncfile]
p_cat = getvar(wrflist, "P", timeidx=ALL_TIMES, method="cat")

print(p_cat)
del p_cat<jupyter_output><empty_output><jupyter_text>2.0.2 Combining via the 'join' method<jupyter_code>p_join = getvar(wrflist, "P", timeidx=ALL_TIMES, method="join")
print(p_join)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>from wrf import ALL_TIMES
p_join = getvar(wrflist, "P", timeidx=ALL_TIMES, method="join", squeeze=False)
print(p_join)
del p_join<jupyter_output><empty_output><jupyter_text>2.0.3 Dictionary Sequences<jupyter_code>wrf_dict = {"label1" : [ncfile, ncfile],
           "label2" : [ncfile, ncfile]}
p_dict = getvar(wrf_dict, "P", timeidx=ALL_TIMES)
print(p_dict)
del p_dict<jupyter_output><empty_output><jupyter_text>2.0.4 Generator Sequences<jupyter_code>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<jupyter_output><empty_output><jupyter_text>2.0.5 Custom Iterable Classes<jupyter_code>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(obj_gen, "P", method="join", squeeze=False)
print(p_obj_gen)

del p_obj_gen
        
import tempfile
import glob
import shutil
import os
from netCDF4 import Dataset
from wrf import getvar, ll_to_xy, CoordPair, GeoBounds


class FileReduce(object):
    def __init__(self, filenames, geobounds, tempdir=None, delete=True, reuse=True):
        """An iterable object for cutting out geographic domains.
        
        Args:
        
            filenames (sequence): A sequence of full paths to the WRF files
            
            geobounds (GeoBounds): A GeoBounds object defining the region of interest
            
            tempdir (str): The location to store the temporary cropped data files. If None, tempfile.mkdtemp is used.
            
            delete (bool): Set to True to delete the cropped files when FileReduce is garbage collected.
            
            reuse (bool): Set to True when you want to resuse the files that were previously converted. *tempdir* 
                must be set to a specific directory that contains the converted files.
        
        """
        self._filenames = filenames
        self._i = 0
        self._geobounds = geobounds
        self._delete = delete
        self._cache = set()
        self._own_data = True
        self._reuse = reuse
        
        if tempdir is not None:
            if not os.path.exists(tempdir):
                os.makedirs(tempdir)
            self._tempdir = tempdir
            if self._reuse:
                self._cache = set((os.path.join(self._tempdir, name) 
                                   for name in os.listdir(self._tempdir)))
        else:
            self._tempdir = tempfile.mkdtemp()

        print ("temporary directory is: {}".format(self._tempdir))
        self._prev = None
        self._set_extents()
    
    def _set_extents(self):
        fname = list(self._filenames)[0]
        with Dataset(fname) as ncfile:
            lons = [self._geobounds.bottom_left.lon, self._geobounds.top_right.lon]
            lats = [self._geobounds.bottom_left.lat, self._geobounds.top_right.lat]
            orig_west_east = len(ncfile.dimensions["west_east"])
            orig_south_north = len(ncfile.dimensions["south_north"])
            
            
            # Note: Not handling the moving nest here
            # Extra points included around the boundaries to ensure domain is fully included
            x_y = ll_to_xy(ncfile, lats, lons, meta=False)
            self._start_x = 0 if x_y[0,0] == 0 else x_y[0,0] - 1
            self._end_x = orig_west_east - 1 if x_y[0,1] >= orig_west_east - 1 else x_y[0,1] + 1
            self._start_y = 0 if x_y[1,0] == 0 else x_y[1,0] - 1
            self._end_y = orig_south_north if x_y[1,1] >= orig_south_north - 1 else x_y[1,1] + 1
            
            self._west_east = self._end_x - self._start_x + 1
            self._west_east_stag = self._west_east + 1
            self._south_north = self._end_y - self._start_y + 1
            self._south_north_stag = self._south_north + 1
            
        
    def __iter__(self):
        return self
    
    def __copy__(self):
        cp = type(self).__new__(self.__class__)
        cp.__dict__.update(self.__dict__)
        cp._own_data = False
        cp._delete = False
        
        return cp
    
    def __del__(self):
        if self._delete:
            shutil.rmtree(self._tempdir)
    
    def reduce(self, fname):
        outfilename = os.path.join(self._tempdir, "reduced_" + os.path.basename(fname))
        
        # WRF-Python can iterate over sequences several times during a 'getvar', so a cache is used to 
        if outfilename in self._cache:
            return Dataset(outfilename)
        
        # New dimension sizes
        dim_d = {"west_east" : self._west_east,
                 "west_east_stag" : self._west_east_stag,
                 "south_north" : self._south_north,
                 "south_north_stag" : self._south_north_stag
                }
        
        # Data slice sizes for the 2D dimensions
        slice_d = {"west_east" : slice(self._start_x, self._end_x + 1),
                   "west_east_stag" : slice(self._start_x, self._end_x + 2),
                   "south_north" : slice(self._start_y, self._end_y + 1),
                   "south_north_stag" : slice(self._start_y, self._end_y + 2)
                  }
        
        with Dataset(fname) as infile, Dataset(outfilename, mode="w") as outfile:
            print ("reduce getting called!")
            
            # Copy the global attributes
            outfile.setncatts(infile.__dict__)

            # Copy Dimensions, limiting south_north and west_east to desired domain
            for name, dimension in infile.dimensions.items():
                dimsize = dim_d.get(name, len(dimension))
                outfile.createDimension(name, dimsize)

            # Copy Variables  
            for name, variable in infile.variables.iteritems():
                
                new_slices = tuple((slice_d.get(dimname, slice(None)) for dimname in variable.dimensions))

                outvar = outfile.createVariable(name, variable.datatype, variable.dimensions)

                outvar[:] = variable[new_slices]

                outvar.setncatts(variable.__dict__)
                
        
        result = Dataset(outfilename)
            
        self._cache.add(outfilename)
            
        return result
            
    
    def next(self):
        if self._i >= len(self._filenames):
            if self._prev is not None:
                self._prev.close()
            raise StopIteration
        else:
            fname = self._filenames[self._i]
            reduced_file = self.reduce(fname)
            if self._prev is not None:
                self._prev.close()
            self._prev = reduced_file
            
            self._i += 1
            
            return reduced_file
    
    # Python 3
    def __next__(self):
        return self.next()

ll = CoordPair(lat=24.0, lon=-87.)
ur = CoordPair(lat=27.0, lon=-84)
bounds = GeoBounds(ll, ur)
reduced_files = FileReduce(glob.glob("/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02*"),
                    bounds, tempdir="/Users/ladwig/mytemp", delete=False, reuse=True)

slp = getvar(reduced_files, "slp")
print(slp)


del (reduced_files)<jupyter_output><empty_output><jupyter_text>3.0 WRF Variable Computational Routines<jupyter_code>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", "wspd_wdir", "wspd_wdir10",
            "uvmet_wspd_wdir", "uvmet10_wspd_wdir"]
#wrf_vars = ["slp"]

vard = {varname: getvar(ncfile, varname, method="cat", squeeze=True) for varname in wrf_vars}
for varname in wrf_vars:
    print(vard[varname])
    print ("\n")<jupyter_output><empty_output><jupyter_text>(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 'to_np' method from wrf.)<jupyter_code>from wrf import to_np
masked_ndarray = to_np(vard["slp"])
print(type(masked_ndarray))
del masked_ndarray
keys = [x for x in vard.keys()]
for key in keys:
    del vard[key]<jupyter_output><empty_output><jupyter_text>3.1 Interpolation Routines  3.1.1 Horizontal Level Interpolation<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>3.1.2 Vertical Cross Section Interpolation<jupyter_code># Pressure using pivot and angle
from wrf import getvar, vertcross, CoordPair

z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
pivot_point = CoordPair((z.shape[-1]-1) // 2, (z.shape[-2] - 1) // 2) 
angle = 90.0

p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle, latlon=True)
print(p_vert)
print ("\n")
del p_vert

# Pressure using start_point and end_point
start_point = CoordPair(0, (z.shape[-2]-1) // 2)
end_point = CoordPair(-1, (z.shape[-2]-1) // 2)

p_vert = vertcross(p, z, start_point=start_point, end_point=end_point, latlon=True)
print(p_vert)
del p_vert, p, z
# Pressure using pivot and angle
from wrf import getvar, vertcross, CoordPair, xy_to_ll

z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
lats = getvar(ncfile, "lat")
lons = getvar(ncfile, "lon")

#print ((lats.shape[-2]-1) / 2)
#print ((lats.shape[-1]-1) / 2)

#print (to_np(lats[529, 899]))
#print (to_np(lons[529, 899]))

#print (to_np(lats[529, 0]))
#print (to_np(lons[529, 0]))

#print (to_np(lats[529, -1]))
#print (to_np(lons[529, -1]))

pivot_point = CoordPair(lat=38.5, lon=-97.5) 
angle = 90.0

p_vert = vertcross(p, z, wrfin=ncfile, pivot_point=pivot_point, angle=angle, latlon=True)
print (p_vert)
print ("\n")

start_lat = lats[(lats.shape[-2]-1)//2, 0]
end_lat = lats[(lats.shape[-2]-1)//2, -1]
start_lon = lons[(lats.shape[-2]-1)//2, 0]
end_lon = lons[(lats.shape[-2]-1)//2, -1]

print (start_lat)
print (end_lat)
print (start_lon)
print (end_lon)

# Pressure using start_point and end_point
start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)

p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)
print(p_vert)
# Pressure using pivot and angle
from wrf import getvar, vertcross, CoordPair, xy_to_ll

z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
lats = getvar(ncfile, "lat")
lons = getvar(ncfile, "lon")

#print ((lats.shape[-2]-1) / 2)
#print ((lats.shape[-1]-1) / 2)

#print (to_np(lats[529, 899]))
#print (to_np(lons[529, 899]))

#print (to_np(lats[529, 0]))
#print (to_np(lons[529, 0]))

#print (to_np(lats[529, -1]))
#print (to_np(lons[529, -1]))

pivot_point = CoordPair(lat=38.5, lon=-97.5) 
angle = 90.0

p_vert = vertcross(p, z, wrfin=ncfile, pivot_point=pivot_point, angle=angle, latlon=True)
print (p_vert)
print ("\n")

start_lat = lats[(lats.shape[-2]-1)//2, 0]
end_lat = lats[(lats.shape[-2]-1)//2, -1]
start_lon = lons[(lats.shape[-2]-1)//2, 0]
end_lon = lons[(lats.shape[-2]-1)//2, -1]

print (start_lat)
print (end_lat)
print (start_lon)
print (end_lon)

# Pressure using start_point and end_point
start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)

levels = [1000., 2000., 3000.]

p_vert = vertcross(p, z, wrfin=ncfile, levels=levels, start_point=start_point, end_point=end_point, latlon=True)
print(p_vert)<jupyter_output><empty_output><jupyter_text>3.1.3 Interpolate 2D Variable to a Line<jupyter_code># T2 using pivot and angle
from wrf import interpline, getvar, CoordPair

t2 = getvar(ncfile, "T2")
pivot_point = CoordPair((t2.shape[-1]-1)//2, (t2.shape[-2]-1)//2) 
angle = 0.0

t2_line = interpline(t2, pivot_point=pivot_point, angle=angle, latlon=True)
print(t2_line, "\n")

del t2_line

# T2 using start_point and end_point
start_point = CoordPair((t2.shape[-1]-1)//2, 0)
end_point = CoordPair((t2.shape[-1]-1)//2, -1)

t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True)
print(t2_line, "\n")

del t2_line

t2 = getvar(ncfile, "T2")
lats = getvar(ncfile, "lat")
lons = getvar(ncfile, "lon")

start_lat = lats[0, (lats.shape[-1]-1)//2]
end_lat = lats[-1, (lats.shape[-1]-1)//2]
start_lon = lons[0, (lons.shape[-1]-1)//2]
end_lon = lons[-1, (lons.shape[-1]-1)//2]

start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)

t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)
print (t2_line)

del t2_line, t2<jupyter_output><empty_output><jupyter_text>3.1.4 Vertical Coordinate Interpolation<jupyter_code>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<jupyter_output><empty_output><jupyter_text>3.2 Lat/Lon to X/Y Routines<jupyter_code>from wrf import xy_to_ll, ll_to_xy 

a = xy_to_ll(ncfile, 400, 200)
a1 = ll_to_xy(ncfile, a[0], a[1])

#print(a)
#print("\n")
#print(a1)
#print("\n")

a = xy_to_ll(ncfile, [400,105], [200,205])
a1 = ll_to_xy(ncfile, a[0,:], a[1,:])

b = ll_to_xy(ncfile, 45.5, -110.8, as_int=True)

# Note: Lists/Dictionaries of files will add a new dimension ('domain') only if the domain is moving
c = xy_to_ll([ncfile, ncfile, ncfile], [400,105], [200,205])
d = xy_to_ll({"label1" : [ncfile, ncfile],
              "label2" : [ncfile, ncfile]}, 
             [400,105], [200,205])

print(a)
print("\n")
print(a1)
print("\n")
print(b)
print("\n")
print(c)
print("\n")
print(d)<jupyter_output><empty_output><jupyter_text>4.0 Plotting with Cartopy<jupyter_code>%matplotlib inline
# 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 to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords

slp = getvar(ncfile, "slp")
smooth_slp = smooth2d(slp, 3)
lats, lons = latlon_coords(slp)

cart_proj = get_cartopy(slp)

fig = plt.figure(figsize=(10,10))
ax = plt.axes(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(), to_np(lons), to_np(lats))
x = xform_coords[:,:,0]
y = xform_coords[:,:,1]

plt.contour(x, y, to_np(smooth_slp), 10, colors="black")
plt.contourf(x, y, to_np(smooth_slp), 10)
plt.colorbar(ax=ax, shrink=.47)

ax.set_xlim(cartopy_xlim(slp))
ax.set_ylim(cartopy_ylim(slp))
ax.gridlines()

# SLP
from __future__ import (absolute_import, division, print_function, unicode_literals)
    
from netCDF4 import Dataset   
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature

from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords, geo_bounds

ncfile = Dataset("/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00")

# Get the sea level pressure
slp = getvar(ncfile, "slp")

# Smooth the sea level pressure since it tends to be noisey near the mountains
smooth_slp = smooth2d(slp, 3)

# Get the numpy array from the XLAT and XLONG coordinates
lats, lons = latlon_coords(slp, as_np=True)

# The cartopy() method returns a cartopy.crs projection object
cart_proj = get_cartopy(slp)

print (cart_proj)

print (get_cartopy(wrfin=ncfile))

bounds = geo_bounds(slp)

print (bounds)

subset = slp[150:250, 150:250]

subset_bounds = geo_bounds(subset)

print (subset_bounds)

file_bounds = geo_bounds(wrfin=ncfile)

print (file_bounds)

# Create a figure that's 10x10
fig = plt.figure(figsize=(10,10))
# Get the GeoAxes set to the projection used by WRF
ax = plt.axes(projection=cart_proj)

# Download and add the states and coastlines
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)

# Make the contour outlines and filled contours for the smoothed sea level pressure.
# The transform keyword indicates that the lats and lons arrays are lat/lon coordinates and tells 
# cartopy to transform the points in to grid space.
plt.contour(lons, lats, to_np(smooth_slp), 10, colors="black", transform=crs.PlateCarree())
plt.contourf(lons, lats, to_np(smooth_slp), 10, transform=crs.PlateCarree())

# Add a color bar
plt.colorbar(ax=ax, shrink=.47)

# Set the map limits
ax.set_xlim(cartopy_xlim(slp))
ax.set_ylim(cartopy_ylim(slp))

# Add the gridlines
ax.gridlines()
# 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, to_np, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords


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)

lats, lons = latlon_coords(ht_500)

cart_proj = get_cartopy(slp)

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(), to_np(lons), to_np(lats))
x = xform_coords[:,:,0]
y = xform_coords[:,:,1]

plt.contour(x, y, to_np(ht_500), 20, cmap=get_cmap("plasma"))
plt.barbs(x[::50,::50], y[::50,::50], to_np(u_500[::50, ::50]), to_np(v_500[::50, ::50]))
plt.colorbar(ax=ax, shrink=.7)

ax.set_xlim(cartopy_xlim(slp))
ax.set_ylim(cartopy_ylim(slp))
ax.gridlines()


# 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, to_np, CoordPair

p = getvar(ncfile, "pressure")
z = getvar(ncfile, "z", units="dm")

pivot_point = CoordPair(z.shape[-1] // 2, z.shape[-2] // 2) 
angle = 90.0

p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle)#, levels=[1000,850,500])

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"))<jupyter_output><empty_output><jupyter_text>Multi-time Moving Domain Files<jupyter_code>import os
from wrf import getvar, ALL_TIMES
from netCDF4 import Dataset as nc

dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest"
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])

p = getvar(ncfiles, "ctt", timeidx=ALL_TIMES)
print (p)
#print (p.attrs["projection"].shape)

print (p.attrs["projection"])

ncfiles[2].variables["XTIME"][:]

p = getvar(ncfiles, "P", timeidx=None, method="cat", meta=True, squeeze=True)

print (p)
print (type(p.coords["Time"]))
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))
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", "cfrac", "uvmet_wspd_wdir",
            "uvmet10_wspd_wdir", "wspd_wdir", "wspd_wdir10"]
#wrf_vars = ["cape_2d"]

vard = {}
for varname in wrf_vars:
    print (varname)
    vard[varname] = getvar(ncfiles, varname, timeidx=None, method="join", squeeze=False)
    
#vard = {varname: getvar(ncfiles, varname, method="join", squeeze=False) for varname in wrf_vars}
for varname in wrf_vars:
    print(vard[varname])
    
# NOTE: Warnings below are due to "join" and the fill values used since the last file only contains 1 time step.
    
import os
from wrf import getvar
from netCDF4 import Dataset as nc

dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest"
ncfilenames = [os.path.join(dir, x) for x in os.listdir(dir) if x.find("_d02_") > 0]
ncfiles = [nc(x) for x in ncfilenames]

# Pressure using pivot and angle
from wrf import getvar, vertcross, CoordPair

timeidx = 0
z = getvar(ncfiles, "z", timeidx, method="join")
p = getvar(ncfiles, "pressure", timeidx, method="join")
pivot_point = CoordPair(z.shape[-1] / 2, z.shape[-2] / 2) 
angle = 40.0

p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle)
print(p_vert)
print ("\n")
del p_vert

# Pressure using start_point and end_point
start_point = CoordPair(0, z.shape[-2]/2)
end_point = CoordPair(-1, z.shape[-2]/2)

p_vert = vertcross(p, z, start_point=start_point, end_point=end_point)
print(p_vert)
del p_vert, p, z
import os
from wrf import getvar
from netCDF4 import Dataset as nc

dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest"
ncfilenames = [os.path.join(dir, x) for x in os.listdir(dir) if x.find("_d02_") > 0]
ncfiles = [nc(x) for x in ncfilenames]

timeidx = None

# T2 using pivot and angle
from wrf import interpline, getvar, to_np, CoordPair

t2 = getvar(ncfiles, "T2", timeidx)
pivot_point = CoordPair(t2.shape[-2] / 2, t2.shape[-1] / 2) 
angle = 90.0

t2_line = interpline(t2, pivot_point=pivot_point, angle=angle, latlon=True)
print(t2_line)
print("\n")


del t2_line

# T2 using start_point and end_point
start_point = CoordPair(t2.shape[-2]/2, 0)
end_point = CoordPair(t2.shape[-2]/2, -1)

t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True)
print(t2_line)
print("\n")

del t2_line, t2
from wrf import getvar

from wrf import xy_to_ll, ll_to_xy 

a = xy_to_ll(ncfiles, 400, 200)
a = xy_to_ll(ncfiles, [400,105], [200,205])
b = ll_to_xy(ncfiles, 45.5, -110.8, as_int=True)

# Note: Lists/Dictionaries of files will add a new dimension ('domain') only if the domain is moving
c = xy_to_ll(ncfiles, [400,105], [200,205])
d = xy_to_ll({"label1" : ncfiles,
              "label2" : ncfiles}, 
             [400,105], [200,205])

print(a)
print("\n")
print(b)
print("\n")
print(c)
print("\n")
print(d)

from glob import glob
from wrf import getvar, ALL_TIMES, geo_bounds, get_cartopy, get_basemap, get_pyngl, cartopy_xlim, cartopy_ylim
from netCDF4 import Dataset as nc

wrf_filenames = glob("/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_*")
ncfiles = [nc(x) for x in wrf_filenames]

slp = getvar(ncfiles, "slp", timeidx=ALL_TIMES)

bounds = geo_bounds(slp)
print (bounds)
print ()

cart_proj = get_cartopy(slp)
print (cart_proj)
print ("\n")

xlims = cartopy_xlim(slp)
print (xlims)
print ("\n")

ylims = cartopy_ylim(slp)
print (ylims)
print ("\n")


bm = get_basemap(slp)
print (bm)
print ("\n")

pyngl_res = get_pyngl(slp)
print (pyngl_res)
print ("\n")<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup