PEP 8 Coding Standards for Python

For standardized python scripting, study this resource...
https://www.python.org/dev/peps/pep-0008/

Adding fields to a structured numpy array created from genfromtxt

Is it possible to add a new field in a numpy.genfromtxt output? Yes it is. I asked this question on stackoverflow and got an answer that works...
http://stackoverflow.com/questions/40182466/is-it-possible-to-add-a-new-field-in-a-numpy-genfromtxt-output/40183748#40183748

The trick is using the numpy.lib.recfunctions library. Use the append_fields function to add a new field to a data structure.

Sorting multiple vectors based on one vector sort order

That title was a mouthful. First, I want to sort an array vector. Then I want to sort other vectors with that same sort order.

I have a vector that contain heights, and with that two other vectors with the corresponding temperature and dew point. I want to sort the height array, and sort the temperature and dew point array in the same manner.

height = np.array([10, 50, 20, 30, 60, 20, 80])
temp = np.array([1, 4, 2, 7, 5, 3, 9])
dwpt = np.array([6, 3, 8, 5, 7, 4, 6])

First, use np.argsort to get the sort order.

sort_order = np.argsort(height)

[out] array([0, 2, 5, 3, 1, 4, 6], dtype=int64)

Now make new arrays of the sorted data with that sort order

sorted_height = np.array(height)[sort_order]
[out] array([10, 20, 20, 30, 50, 60, 80])

sorted_temp = np.array(temp)[sort_order]
[out] array([1, 2, 3, 7, 4, 5, 9])

sorted_dwpt = np.array(dwpt)[sort_order]
[out] array([6, 8, 4, 5, 3, 7, 6])

Thus, you see that we have sorted the height array from lowest to highest, and the temperature and dew point arrays have been sorted, preserving their order with the corresponding height.

Spanish Fork Climatology

I just used the MesoWest API to create a climatology of weather variables at my weather station in Spanish Fork. This shows the max, mean, and min values for various variables of each month of the year. The data collected ranges from  2013 to 2016.

See code: https://github.com/blaylockbk/Ute_MesoWest/blob/master/station_month_climatology.py

The freezing temperature is marked in a dashed light blue.

And looking at the daily averages (again this is only a few years of data, so in no way an official "climatological" average of temperature on each day)

See code: https://github.com/blaylockbk/Ute_MesoWest/blob/master/station_day_climatology.py

We haven't had a day below freezing after the first few days of May.

MetPy pythons gift to meteorologists :)

HRRR 700mb wind field and vorticity over Utah

I just installed MetPy and started using it to calculate and plot maps of HRRR vorticity! This will be a very useful package in the future :)

Python Transparent Colormap

Transparent color maps are very handy when plotting data on top of other data. If found this little snipped of code...

import matplotlib.colors as mcolors

colors = [(1,0,0,c) for c in np.linspace(0,1,100)]
cmapred = mcolors.LinearSegmentedColormap.from_list('mycmap', colors, N=5)
colors = [(0,0,1,c) for c in np.linspace(0,1,100)]

cmapblue = mcolors.LinearSegmentedColormap.from_list('mycmap', colors, N=5)

Then you can assign the colormap in your plot

plt.contourf(x,y,z,cmap=cmapred)

plt.colormesh(x,y,z,cmap=cmapblue)

You may need to adjust the levels to get the scale right in the contours.

Here's an example of my transparent colormap showing pollution tracers blowing around northern Utah from Salt Lake County...

Below is how I made the above 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
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
#import mayavi.mlab as mlab #Used for 3d Plotting
import matplotlib.colors as mcolors
from functions_domains_models import *

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

# list of files
#d02 = ['WTDE8Z~0','W2TFAX~5','W1PEFC~6','WKBXUM~Z','WAMGEL~S','WWZHSX~P']
#d02 = ['WKBXUM~Z']
d02 = ['auxhist24_d02_2015-06-17_00:00:00']
spat= 'auxhist23_d02_2015-06-17_00:00:00'
for output_file in d02:
# Open file in a netCDF reader
directory = '/uufs/chpc.utah.edu/common/home/horel-group4/model/bblaylock/WRF3.7_lake303_ember/WRFV3/test/em_real/'
out_dir = '/uufs/chpc.utah.edu/common/home/u0553130/public_html/MS/tracercomposite/'
if not os.path.exists(out_dir):
os.makedirs(out_dir)
wrf_file_name = output_file #NC files are named differently in a windows system
print 'opening', directory+wrf_file_name
nc = netcdf.netcdf_file(directory+wrf_file_name,'r')
nc_spatial = netcdf.netcdf_file(directory+spat,'r')



#-------------------------------------------
# Creat the basemap only once per domain. (This significantly speeds up the plotting speed)
#-------------------------------------------

# 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

""" #Full map
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)
"""
domain = get_domain('salt_lake_valley')
top_right_lat = domain['top_right_lat']+1
top_right_lon = domain['top_right_lon']+1
bot_left_lat = domain['bot_left_lat']-1
bot_left_lon = domain['bot_left_lon']-1

## Map in cylindrical projection (data points may apear skewed)
m = Basemap(resolution='i',projection='cyl',\
llcrnrlon=bot_left_lon,llcrnrlat=bot_left_lat,\
urcrnrlon=top_right_lon,urcrnrlat=top_right_lat,)


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

# Define the colormaps
colors = [(1,0,0,c) for c in np.linspace(0,1,100)]
cmapred = mcolors.LinearSegmentedColormap.from_list('mycmap', colors, N=8)
colors = [(0,0,1,c) for c in np.linspace(0,1,100)]
cmapblue = mcolors.LinearSegmentedColormap.from_list('mycmap', colors, N=8)


for t in np.arange(0,int(np.shape(nc.variables['Times'])[0])):
plt.cla()
plt.clf()
plt.close()

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

time = ''.join(nc.variables['Times'][t])
time_dt = datetime.strptime(time,'%Y-%m-%d_%H:%M:%S')
time_str= datetime.strftime(time_dt,'%Y-%m-%d %H:%M:%S UTC')
time_str_local=datetime.strftime(time_dt-timedelta(hours=6),'%Y-%m-%d %H:%M:%S MDT')
time_file=datetime.strftime(time_dt,'%Y%m%d%H%M%S')
HGT = nc_spatial.variables['HGT'][0,:,:] #topography
landmask = nc_spatial.variables['LANDMASK'][0,:,:]
plumeN = nc.variables['N_SLV'][t][0]
plumeN_composite = np.sum(nc.variables['N_SLV'][t], axis=0)
plumeS = nc.variables['S_SLV'][t][0]
plumeS_composite = np.sum(nc.variables['S_SLV'][t], axis=0)

#d = nc.variables['tr17_2'][0]
#mlab.contour3d(d)


# 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(XLONG[0,::thin,::thin],XLAT[0,::thin,::thin])

# Set universal figure margins
width = 10
height = 8
plt.figure(t+1)
print 'Plotting',time_str
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)

######################################################

plt.contour(x,y,landmask,levels=[1])
plt.contourf(x,y,HGT,cmap=plt.get_cmap('binary'))
#plt.pcolormesh(x,y,plume,cmap='PuOr', vmin=0,vmax=.01)
plt.pcolormesh(x,y,plumeN_composite,cmap=cmapred,vmax=8,vmin=0)
plt.pcolormesh(x,y,plumeS_composite,cmap=cmapblue,vmax=8,vmin=0)
cbar_loc = plt.colorbar(shrink=.8,pad=.01,ticks= np.arange(0,8,1))
cbar_loc.ax.set_ylabel('Tracers [units of stuff]')
plt.contour(x,y,landmask, [0,1], linewidths=1, colors="k")
#plt.contour(x,y,HGT)

xs,ys = m(-111.97,40.78) #buffr sounding coordinates
plt.scatter(xs,ys, s=70, c='w')
plt.title(time_str+'\n'+time_str_local, bbox=dict(facecolor='white', alpha=0.65),\
x=0.5,y=.90,weight = 'demibold',style='oblique', \
stretch='normal', family='sans-serif')
plt.savefig(out_dir+'TracersComposite_'+time_file+'.png',bbox_inches="tight")
print 'Saved',time_str

#plt.show()

python colormaps

When creating plots it is very important to consider the kind of colormap you wish to display the data in. Some colormaps will emphasize certain features more than others. But other color schemes may be deceiving, particularly the rainbow color scheme. See an example here.

Utah terrain with
viridis color map in matplotlib 2.0

Note: In Matplolib Version 2 the default colormap is a green shade called 'viridis' which is much better than jet (see here).

Python's matplotlib module has many preloaded colormaps you can use in your figures. You can look at the available colormaps here: http://matplotlib.org/examples/color/colormaps_reference.html, but I find that looking at just the color bars makes it difficult to imagine what a real plot might look like. So...I created a terrain plot for Utah at 3 km resolution with each python colormap to help visualize what each of the plot would look like.

To set your plot with the desired colormap, simply use the argument:
cmap=plt.get_cmap('name_of_colormap')

For example:
plt.pcolormesh(x,y,height,cmap=plt.get_cmap('Accent'))

Below are all the available colormaps in matplotlib. It is possible to create your own more info on creating your own here.

Matplotlib Default Colormaps

 

Custom Terrain Colormap
It is possible to create your own colormap (click  here).
A custom terrain colormap created using this custom colormap: here.

Luke's Colormaps

From github https://github.com/lmadaus/old_wrf_plotting_scripts/blob/master/coltbls.py. You'll load these differently than the default colormaps, but the following are the different options you can use...