Quick WRF Plot

Plotting WRF Data in python without NetCDF4:

Scrapped this code from places on line. Reads in netcdf data with scipy. Plots specific humidity with the lake mask outlined on top. See blog post here for details on the plot.

# Brian Blaylock
#
# Plotting WRF netCDF output
# Parts and Pieces of this code is from Luke Madaus (University of Washington)


import sys,getopt
#from netCDF4 import Dataset # use scipy instead
from scipy.io import netcdf
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import os
from datetime import datetime, timedelta
import coltbls as coltbls
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.axes_grid import make_axes_locatable
import matplotlib.axes as maxes


# Running Script example:
# $ python my_WRF_map.py wrfout_d0x_YYYY-MM-DD_HH:MM:SS

# Open file in a netCDF reader
directory = ''
wrf_file_name = 'M11KS3~E.NC' #NC files are named differently in a windows system
nc = netcdf.netcdf_file(wrf_file_name,'r')

# Grab these variables for now
landmask = nc.variables['LANDMASK'] # HRRR metgrid landmask




#-------------------------------------------
# BEGIN ACTUAL PROCESSING HERE
#-------------------------------------------

# x_dim and y_dim are the x and y dimensions of the model
# domain in gridpoints
x_dim = nc.dimensions['west_east']
y_dim = nc.dimensions['south_north']

# Get the grid spacing
dx = float(nc.DX)
dy = float(nc.DY)

width_meters = dx * (x_dim - 1) #Domain Width
height_meters = dy * (y_dim - 1) #Domain Height

cen_lat = float(nc.CEN_LAT)
cen_lon = float(nc.CEN_LON)
truelat1 = float(nc.TRUELAT1)
truelat2 = float(nc.TRUELAT2)
standlon = float(nc.STAND_LON)

# Draw the base map behind it with the lats and
# lons calculated earlier
m = Basemap(resolution='i',projection='lcc',\
width=width_meters,height=height_meters,\
lat_0=cen_lat,lon_0=cen_lon,lat_1=truelat1,\
lat_2=truelat2)

# This sets the standard grid point structure at full resolution
# Converts WRF lat and long to the maps x an y coordinate
x,y = m(nc.variables['XLONG_M'][0],nc.variables['XLAT_M'][0])


# This sets a thinn-ed out grid point structure for plotting
# wind barbs at the interval specified in "thin"
thin = 5
x_th,y_th = m(nc.variables['XLONG_M'][0,::thin,::thin],\
nc.variables['XLAT_M'][0,::thin,::thin])

# Set universal figure margins
width = 10
height = 8

plt.figure(figsize=(width,height))
plt.rc("figure.subplot", left = .001) #gets rid of white space in plot
plt.rc("figure.subplot", right = .999)
plt.rc("figure.subplot", bottom = .001)
plt.rc("figure.subplot", top = .999)

F = plt.gcf() # Gets the current figure

m.drawstates(color='k', linewidth=1.25)
m.drawcoastlines(color='k')
m.drawcountries(color='k', linewidth=1.25)

height = nc.variables['HGT_M'][0,:,:]
landmask = nc.variables['LANDMASK'][0,:,:]
spechumd = nc.variables['SPECHUMD'][0][0,:,:]

plt.figure(1)
plt.title("HRRR Land Mask and Specific Humidity\nAugust 10, 2015 23:00z")
plt.pcolormesh(x,y,landmask)
plt.pcolormesh(x,y,spechumd*1000, cmap = plt.cm.cool)
cbar_loc = plt.colorbar()
cbar_loc.ax.set_ylabel('Specific Humidity g/kg')
plt.contour(x,y,landmask, [0,1], linewidths=1, colors="k")
#plt.contour(x,y,height)

xs,ys = m(-111.97,40.78) #buffr sounding coordinates
plt.scatter(xs,ys, s=70, c='w')







plt.show()