Unidata Python Training

Some good resources to learn Python

https://unidata.github.io/online-python-training/

Plotting radar data with MetPy, pyproj, Basemap

MetPy radar plots

The MetPy python package is really helpful for atmospheric scientists which allows you to plot radar data (specifically Level 3 data) downloaded from the THREDDS server whose files are in the .nids format. How to use MetPy to plot these data is explained on the MetPy docs, and you'll create a plots of reflectivity and radial velocity like this...

There are two things the example doesn't tell you...

1) The data value units are arbitrary. Notice that the values of dbz color scale is the same as the values in the radial velocity scale on the right. What the MetPy example doesn't tell you is that you need to use the map_data instance to convert the data to appropriate units.

f = Level3File(FILE)

# Pull the data out of the file object (this stuff is in 8-bit)
datadict = f.sym_block[0][0]

# Turn into an array, then mask
data = f.map_data(datadict['data']) # I'm confident this maps the data to floats with units of m/s
data = ma.array(data)
data[np.isnan(data)] = ma.masked

2) You don't get a sense of where the radar and its observations are located in relation to places, just a distance from the center. You can convert the azimuth and range into a latitude and longitude by using the pyproj module. Then you can use the lon and lat variables to plot the data on a basemap. Given the center latitude, center longitude, azimuth (degrees from north), and the range from center lat/lon, you can calculate the latitude and longitude of a distant point like this...

from pyproj import Geod

# Grab azimuths and calculate a range based on number of gates
az = np.array(datadict['start_az'] + [datadict['end_az'][-1]])
rng = np.linspace(0, f.max_range, data.shape[-1] + 1)

# Convert az, range to a lat/lon
g = Geod(ellps='clrk66') # This is the type of ellipse the earth is projected on. 
                         # There are other types of ellipses you can use,
                         # but the differences will be small
center_lat = np.ones([len(az),len(rng)])*f.lat    
center_lon = np.ones([len(az),len(rng)])*f.lon
az2D = np.ones_like(center_lat)*az[:,None]
rng2D = np.ones_like(center_lat)*np.transpose(rng[:,None])*1000
lon,lat,back=g.fwd(center_lon,center_lat,az2D,rng2D)
#...(add code to create map)
m.pcolormesh(lon,lat,data,cmap="BrBG_r",vmax=10,vmin=-10)

You can replace the xlocs and ylocs in the MetPy example with the lon and lat variables, and use those lons/lats to plot on a basemap.
(Note: back_az is the azimuth from the point back to the center location)

TDWR plot on Basemap

I like to look at the TDWR base radial velocity scans to look for lake breezes in the Salt Lake Valley. Below is plotted the radial velocity on June 19, 2015 at 00:24 UTC using MetPy without converting to correct units using map_data and has not been plotted on a basemap (two color schemes used).

both plots are the same, but the color bar is different.

It is so much nicer to view radar data on a basemap plot, as shown below, so you get a sense of there the observations are in space. You can see a lake breeze boundary in the center of the valley (below center of the image where color changes from dark green to light green). I have also used the map_data instance to covert the data to units of ms-1. Each plot is the same with different color schemes. Left is the Brown/BlueGreen that I like, and right is the NWS Red/Green scheme built into MetPy.

 

The full code to do this is on github:
Plotting TDWR Radar data:
https://github.com/blaylockbk/pyBKB_v2/blob/master/TDWR/plot_TDWR_reflectivity.py

Plotting NEXRAD Level 2 data:

https://github.com/blaylockbk/pyPlots_v2/blob/master/NEXRAD_II/plot_NEXRAD_II_data.py

Extra things you can put on the map are a shape file of the roads, and weather observations from the MesoWest API. Let me know if you find this useful :)

Plot a psudo-sounding with MesoWest Data

In the Salt Lake City area there are many weather stations at elevations higher than the valley floor. Farnsworth Peak is such a weather station at almost 700 hPa. I wondered what it would look like to plot station weather data on a skew-t graph, and see how it compares to rawinsonde balloons.

Here I use the MesoWest API and the Wyoming Sounding data to plot these graphs. I use the SkewT 1.1.0 module for Python to create the SkewT graph. The balloon observations are in green and the real observations are plotted in blue with the station id next to the point. A red line is drawn to create the "psudo-sounding" from the temperature observations at the different elevations. It doesn't match up exactly, except, for FWP, which seems to be an accurate representation of the temperatures of the atmosphere at that elevation.

See more of the code on GitHub: https://github.com/blaylockbk/pyBKB and https://github.com/blaylockbk/Ute_other/blob/master/MesoWest_psudo_sounding.py

## Brian Blaylock

## 27 June 2016

# A lot of the stations in Utah are at elevations.

# I want to plot station data on a skew-t graph to get a 

# psudo vertical profile for any time of day.

import matplotlib.pyplot as plt

import numpy as np

from datetime import datetime, timedelta

import sys

sys.path.append('B:\pyBKB')                                        #for running on PC (B:\ is my network drive)

from BB_MesoWest import MesoWest_stations_radius as MW

from skewt_900_700 import SkewT #this version plots between 900-700 hPa

import wyoming_sounding

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

request_time = datetime(2016,6,29,12,0)

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

MWdatetime = request_time.strftime('%Y%m%d%H%M')

a = MW.get_mesowest_radius_winds(MWdatetime,'10')

p = a['PRESS']>0 #pressure is in these indexes

pres = a['PRESS'][p]/100

temp = a['TEMP'][p]

hght = a['ELEVATION'][p]

dwpt = a['TEMP'][p]*np.nan #don't get the dewpoint now

wd   = a['WIND_DIR'][p]

ws   = a['WIND_SPEED'][p]

stns = a['STNID'][p]

# sort all that data based on pressure

from operator import itemgetter

pres,temp,hght,dwpt,wd,ws,stns = zip(*sorted(zip(pres,temp,hght,dwpt,wd,ws,stns),key=itemgetter(0)))

data = {'drct':wd,

        'dwpt':dwpt,        #don't plot dwpt right now

        'hght':hght,

        'pres':pres,

        'sknt':ws,

        'temp':temp}

S=SkewT.Sounding(soundingdata=data)

S.plot_skewt(title='psudo sounding\n'+MWdatetime,color='r',lw=.8,parcel_type=None,

             pmin=700.,pmax=900.)

for i in range(0,len(stns)):

    if temp[i]>0: #wont accept nan values

        print temp[i],pres[i],stns[i]

        S.plot_point(temp[i],pres[i])

        S.plot_text(temp[i],pres[i],str(stns[i]))

        

# Add real sounding to plot

# If request_time is in the afternoon, then request the evening sounding

if request_time.hour > 18:

    # Get the evening sounding (00 UTC) for the next day

    next_day = request_time + timedelta(days=1)

    request_sounding = datetime(next_day.year,next_day.month,next_day.day,0)    

elif request_time.hour < 6:

    # Get the evening sounding (00 UTC) for the same

    request_sounding = datetime(request_time.year,request_time.month,request_time.day,0)    

else:

    # Get the morning sounding (12 UTC) for the same day

    request_sounding = datetime(request_time.year,request_time.month,request_time.day,12)

print "SLC sounding for", request_sounding

a = wyoming_sounding.get_sounding(request_sounding)

balloon = {'drct':a['wdir'],

           #'dwpt':a['dwpt'],

           'hght':a['height'],

           'pres':a['press'],

           'sknt':a['wspd'],

           'temp':a['temp']}

# Add to plot

T = SkewT.Sounding(soundingdata=balloon)

S.soundingdata=T.soundingdata

S.add_profile(color='g',lw=.9,bloc=1.2)

             

WRF Browser for Windows

Again, not a python post, but some good news for WRF modelers who use Windows PCs.


Introducing WRF Browser, a Windows application you can load and view WRF output and netcdf variables on a map or export to CSV. http://www.enviroware.com/WRF_Browser/
A few tips unique to my needs:
1) I like to look at terrain and land use fields. This is especially nice when you can do this on Google Earth. You'll need to create your own colormap, or as they call it "ColorRamp". Follow the directions in the online manual for creating a custom colormap. I've made one called MODIS_21_Landuse.pg and put it in the directory C:\ProgramData\Enviroware\WRF_Browser\ColorRamps\ 

This matches my Python ColorMap for MODIS land use

1         0 102   0 255
2         0 102  51 255
3        51 204  51 255
4        51 204 102 255
5        51 153  51 255
6        76 178   0 255
7       209 104  31 255
8       188 181 105 255
9       255 214   0 255
10        0 255   0 255
11        0 255 255 255
12      255 255   0 255
13      255   0   0 255
14      178 230  77 255
15      255 255 255 255
16      233 233 179 255
17      127 178 255 255
18      219  20  59 255
19      247 127  80 255
20      232 150 122 255
21        0   0 224 255



The KML output looks like this...

This is the same land use scheme used in the HRRR model, except that the Size of the Great Salt Lake has been modified to it's present day lake level.
One thing I notice is that urban sprawl in the southern end of Utah County in Spanish Fork (my hometown), Salem, and Payson is not captured in the MODIS land use categories. This may have some small effects on how well the HRRR produces forecasts in these cities (as far as I can tell, land use and surface characteristics for a city is a minor issue with the HRRR compared to it's other problems).

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.