This page demonstrates Python tips and tricks that I use in my everyday programming as an atmospheric science graduate student.
-Brian Blaylock

Tuesday, August 2, 2022

Grid Spacing vs. Resolution

If you spend enough time with meteorologist, you'll inevitably hear criticism against using the non-interchangeable terms "grid spacing" and "resolution" to describe the spacing between grid points in a numerical weather model. Just last week I reviewed a paper where one of the other reviewers blasted the authors for incorrectly using the terms interchangeably. I cringed at the comment. What drives me nuts about this is that the reviewer didn't even explain why the terms are not interchangeable. The supposed incorrect use of the term model "resolution" in place of "grid spacing" will go on forever because those who are deeply offended by those who interchange the term rarely explain how the two terms are different.  

What is the difference?

In a March 2000 BAMS commentary, Lewis D. Grasso explains: 

The two terms  grid  spacing and resolution refer to two different length scales that characterize a grid configuration. Because of this, they are not interchangeable in linear and nonlinear numerical  models. It is not possible to resolve a wave on the scale of one grid spacing in any spatial direction.  Since 2Δx and 3Δx waves are removed to prevent nonlinear instability, waves on the scale of at least 4Δx may be resolved. Other terms that may be used in place of grid spacing are grid interval, grid length, and grid increment

There you have it; resolution is not the correct term to describe a model's grid spacing because resolution refers to the length scale that waves are resolved in the model.

More resources on the topic


Friday, January 15, 2021

Documentation for my Packages

  I finally got around to consolidating some of my more useful Python code and generated Sphinx documentation.



Please check out the following packages:

  • GOES-2-go
    • Download GOES-16 or GOES-17 data from Amazon Web Services (ABI and GLM data).
    • RGB recipes for ABI images (Natural Color, Air Mass, Water Vapor, etc.)
  • SynopticPy:
    • Retrieve mesonet data from the Synoptic API
  • HRRR-B:
    • Download HRRR data from the HRRR archives at the University of Utah's Pando archive, Google Cloud Platform, and Amazon Web Services.

Friday, November 15, 2019

Pygrib now on Windows via anaconda conda-forge

I was very excited to hear this news that pygrib is now available on Windows...



I gave it a try, and it seems to open HRRR grib2 files just fine. This means my HRRR_Pando scripts will work on a Windows machine now. Happy day!

Wednesday, July 10, 2019

How to download specific variables from HRRR GRIB2 files from the Pando archive

Don't have much time, but wanted to paste in a reply to an email concerning downloading specific variables from the HRRR archive on Pando.

The person wanted to download a few variables from files on the archive, rather than the full file. They also wanted to get data from a specific latitude/longitude location in the model grid.

-------------

Unfortunately, it isn't possible to pluck out specific lat/lon locations from the grid before you download the data. The best thing to do would be to download the full grid for the variables you want.
I have described how to do this on the "Scripting Tips" page under the "cURL and wget" and "Python" tabs.

The approach here is that by knowing the byte range for the grid of a single variable (as found in the .idx file), you can download just that chunk from the file (~1 MB) rather than the full file (>125 MB).
An example .idx file for the model analysis (F00) valid at 0000 UTC on 10 July 2019 is here: https://pando-rgw01.chpc.utah.edu/hrrr/sfc/20190710/hrrr.t00z.wrfsfcf00.grib2.idx
This file corresponds to the full GRIB2 file here: https://pando-rgw01.chpc.utah.edu/hrrr/sfc/20190710/hrrr.t00z.wrfsfcf00.grib2

In this example, if you want surface temperature, you would search the .idx file for the line with the variable TMP:surface
The grid for this variable begins on byte 31889605 and the next variable begins on byte 33367664.
Thus, you can download just that grid, which returns a valid GRIB2 file with one layer, with a curl command:
curl -o 20190710_00zf00_TMPsurface.grib2 --range 31889605-33367664 https://pando-rgw01.chpc.utah.edu/hrrr/sfc/20190710/hrrr.t00z.wrfsfcf00.grib2

Because we don't archive the HRRR on the native grids we don't have the data you need at specific model levels. But we do archive the surface file forecasts (F00-F18) and the pressure file analyses (F00) that interpolate the data on specific levels. These still might be useful for what you want to do.

You can browse the variable names and descriptions here: (though, these pages might not be a complete list)

I listed the variables you might be interested in below:


-geopotential (at surface and at first level of model)

HGT:surface
(no HGT at first model level)

-wind speeds (u and v components at first level)

UGRD:10 m
VGRD:10 m
UGRD:80 m
VGRD:80 m
WIND:surface This is the maximum wind speed at 10 m for the previous hour
GUST:surface This is essentially the maximum wind speed in the lowest XX layers (not sure exactly how gust is calculated)
(again, wind is not available at first model level)
-snow depth
SNOD:surface
SNOWC:surface (snow cover)

-temperature (at surface and at first level)

TMP:surface
TMP:2 m
(Not available at first level)
-solar irradiance
DSWRF:surface (downward shortwave radiation flux)
-relative humidity
RH:2 m
-cloud fraction
TCDC:entire atmosphere (total cloud cover)
-hourly precipitation amount
This is a tricky one because precip is an accumulated quantity, thus there isn't a value for the analysis (F00). Instead, you would have to look at accumulated precipitation over some period of time. Some people like to look at the precip accumulated during the F00-F01 period as the "best guess"
For example:
APCP:surface:0-2 hour acc fcst Is the accumulated precip between F00-F02
APCP:surface:1-2 hour acc fcst Is the accumulated precip between F01-F02

If you target the specific variables you want, you will end up downloading much less data than downloading the full files. This will likely save a lot of download time and bandwidth (on your end and our end). It would probably be best to execute a curl command for each variable separately (and for each hour of the day).

I often have to pluck out values at specific lat/lon locations, too. My method is available on github, but there are probably many other ways to do the same thing.

I hope this helps!

Brian Blaylock
Ph.D. Candidate
Atmospheric Science
University of Utah

Thursday, May 30, 2019

RAP and HRRR model domain boundaries

I made a neat map of the RAP (version 4) and HRRR model domain boundaries with Cartopy. I think it looks nice...


Tuesday, May 21, 2019

Cartopy: add NEXRAD mosaic image to figure

This Notebook shows how I add NEXRAD mosaic images from Iowa State University to my Cartopy maps...

https://github.com/blaylockbk/pyBKB_v3/blob/master/BB_maps/cartopy_NEXRAD-mosaic-from-Iowa.ipynb


Slowly, like a sloth, I'm using cartopy more and more.

Monday, May 20, 2019

Add an ArcGIS Map Service Image to a Cartopy axes.

I am a very slow Cartopy adopter. I learned Basemap years ago and am very comfortable with it and continue to hold tight to what I know. However, as I dabble more with Cartopy, I can see how it is superior.

For the last two years, one of my many reasons for sticking with Basemap despite almsot everyone telling me to stop, that the ease of adding arcGIS images as a map background. In Basemap, it is a simple m.arcgisimage(service='World_Shaded_Relief'). It took me a long time and many attempts to figure out how to do this in Cartopy. It does take a few extra steps, but it is possible, and even without the ease of use, I am now believe Cartopy does this better than Basemap.

In Basemap, you are limited to adding arcGIS images to cylindrical projection, but in Cartopy you can add these images to any projection. Incredible!

See the full notebook here, https://github.com/blaylockbk/pyBKB_v3/blob/master/BB_maps/cartopy_arcgisimage.ipynb.


Additional Resources:

NOTE: When you set the extent with ax.set_extent, the images are fine until you cross the prime meridian, then the images are distorted as if there aren't enough tiles. The only way I've been able to prevent this is to keep the image in the same projection as the tiles, a.k.a. the Mercator projection (see the second stack overflow answer linked above).

import cartopy.io.img_tiles as cimgt

fig=plt.figure(figsize=(10,5))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.coastlines()

url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}.jpg'

image = cimgt.GoogleTiles(url=url)

ax.add_image(image, 1)




Wednesday, April 17, 2019

GOES-16 and GOES-17 Overlay

I have successfully plotted GOES-16 and GOES-17 true color images together on the same plot. GOES16 was plotted first, then the GOES17 layer was added on top.

Details on GitHub:
https://github.com/blaylockbk/pyBKB_v3/blob/master/BB_GOES/Both_GOES16_and_GOES17.ipynb


Monday, April 23, 2018

Tips for Publication-worthy Matplotlib Figures

Matplotlib is the graphical plotting package for Python and is heavily used to create. Here are a few tips to create a publication-quality image in Python. There is nothing worse than seeing a poorly made, blurry, or crowded figure when you read a journal article.

The American Meteorological Society has a document that explains best practices to make your figures. These are good rules of thumb to follow.

  1. Create each figure in its own Jupyter Notebook. I recently learned why Jupyter Notebooks are so amazing...it's a great way to document and cleanly see each step. You can easily share notebooks and upload them to GitHub. They also give you the power to make minor changes and do figure 'touch-ups'. These touch-ups are what will make your figures look amazing. 
  2. Pay attention to the dpi: dpi is the 'dots per inch'. A high dpi will make your images look crisp. Generally, anything less than 150 is ideal for viewing on a web, but for print, you want at least 300 dpi. If you figure is just black and white, use at least 1000 dpi. If it includes shades of grey or colors it is acceptable to bump the dpi down to 600. 
  3. Do not manually resize your picture after it has been created. For example, do not make a large picture, load it into PowerPoint, and then shrink the image. This will distort the image, change the dpi, make lines too thin, and cause all sorts of asthetic problems. It is best to make the image the size you want it to be from the beginning.
  4. Pay careful consideration to label, titles, and legends. You want these to be readable, so pay attention to their font size. In fact, most of the time you can get away with not using a figure title--that is the purpose of the figure caption.
  5. Do not use the `Jet` color map. Jet is just horrible. If you need some convincing, then watch this. You can also read a million other articles why Jet is horrible.
  6. Use colors that make sense. If you are showing two diverging values, like negative-ness and positive-ness of a value, use a diverging color map. If you are illustrating hot temperatures, use red.
  7. Consider the right marker shape.
In general, make your figures look good. It is worth spending extra time to make your figures look good. After all, the pretty pictures in your article are the most-looked-at part of your paper.

Tuesday, November 7, 2017

More Python Color maps: cmocean

Additional, built-in colormaps can be found by importing cmocean
http://matplotlib.org/cmocean/

import cmocean
plt.plot(data, cmpa=cmocean.cm.temperature)

Options:
  • temperature
  • salinity
  • par
  • gray
  • oxygen
  • bathymetry
  • density
  • chlorophyll
  • CDOM
  • turbidity
  • speed
  • waveheight
  • waveperiod
  • phase
  • freesurface
  • velocity
  • vorticity
Don't use the Jet colormap!
http://www.research.ibm.com/people/l/lloydt/color/color.HTM

Monday, November 6, 2017

python enumerate()

I have seen the enumerate() function many times in other peoples scripts, but I have never known what it does exactly, and why someone would prefer using it.

Enumerate means to "mention a list of items, one by one." So, it would be most useful in loops over a list.

The built in enumerate() function lists the items in a list along side their index number:

order_list = ['fist', 'second', 'third']
for i in enumerate(order_list):
   print i

(0, 'fist')
(1, 'second')
(2, 'third')


This is useful for getting the index values of a list:

for i, order in enumerate(order_list):
    print i, order

0 fist
1 second
2 third


Thus, this function is a better way of looping over a certain list:

for i in range(len(order_list)):
    print i, order_list[i]

0 fist
1 second
2 third


Why is enumerate() better than range(len())? Because it is one function instead of two nested functions!

You may also iterate through each item of a list without a loop by using the next function:

items = enumerat(order_list)
items.next()

items.next()

items.next()


See more here: http://www.juniordevelopercentral.com/python-enumerate-function/

Friday, August 4, 2017

Lessons Learned using the Open Science Grid

I have been playing with the Open Science Grid, a high throughput computing system. I am using the system to compute statistics of HRRR data over the last 2.5 years stored in the CHPC Pando archive system. The archive is currently 35TB, and growing every day.

OSG Limits: 102,400 MB and 500,000 Files

The OSG appears to run faster in the early morning. Maybe less people are on it.

Downside is the file transfer time if you need to transfer lots of files to your computer. I processed three years of HRRR data, generating almost 700 GB of statistical data for a single variable in 4 hours, but to transfer all those 8,700+ files to my home institution took 1 day and 7 hours. Still, that is faster than running the job on our local nodes, which took 7.5 days to complete the same job on wx4 which has 8 cores.

My third run using the Open Science Grid I calculated 30-day running percentiles for every hour of the year of composite reflectivity, REFC:entire atmosphere, from the HRRR data. It took 2 hours and 10 minutes, running 8,784 jobs (one job for every hour of the year including leap year). Some of the jobs failed (about 10-15) for an unknown reason. I simply dealt with it by re-running the job locally and that didn't take too long. I created 331 GB of data (~43 MB per file), a little smaller than the temperature or wind grids.

I've seen over 455 jobs run simultaneously on the OSG.

File Size

LatLon 18MB
TMP: 2 m  | 684 GB (~82 MB/file) | All Map Images: 8.87 GB | Transfer: 6.5 hours (8 cores) |
WIND:10 m | 704 GB (~87 MB/file) | All Map Images: 8.98 GB | Transfer: 2.5 hours (32 cores) |
DPT:2 m | _ GB (~84 MB/file) | All Map Images: _GB | Transfer: _ hours (32 cores) |
REFC:entire | _ GB (~_ MB/file) | All Map Images: _GB | Transfer: _ hours (32 cores) |

Should have 8,784 jobs

Since I still get errors when running on OSG, I have to do a little babysitting in my workflow:
  1. Create condor_submit script on OSG
    • a python script helps me with this
  2. Submit the job on the OSG:
    • condor_submit job.submit
  3. Copy the files to local CHPC storage using scp
    • Again, another python script loops through all available files
  4. Rerun the statistics calculations locally for files that are missing.
    • For a reason I haven't investigated yet, some files are not created on the OSG. I rerun these locally because running 10-15 files doesn't take as much time as the 8,000+ that are run on OSG.
  5. Rerun bad files
    • Some of the files returned from OSG are incomplete, so I check the file sizes of all the files and when I see one that is smaller than the rest I rerun the job locally.
  6. Remove files from OSG
    • These jobs produce a lot of data, over 500% of my allocation on the OSG, so I quickly remove those files.

A note about saving HDF5 files:
It is best to save each array on "top level" of the HDF5 file. My first iteration of creating these files I stored all my percentile calculations as a 3D array, requiring to get a slab of data to index the individual levels, i.e. file['percentile values'][0][10, 10] to get a point of data. It is much better and more efficient to grab data stored in 2D arrays, on the "top level" like my mean array is stored, i.e. file['mean'][10, 10]. With multiprocessing, I created a time series of the mean data at a point in 30 seconds, but it took 33 minutes to create a time series at a point for the percentile data!!! I'm pretty sure this is due to the added dimension of the percentile variable.

Each HRRR file has a grid size of 1059x1799 (1.9 million grid points), and 136 different meteorological variables. That is a lot of data that needs to be sifted through.

The deal with using the OSG is that your jobs need to be embarrassingly parallel. Thus, I have had framed the research questions I want answered in an embarrassingly parallel framework.
In my specific case, I have chosen to compute hourly statistics for each month--a job that can be done in 288 parts (12 months * 24 hours). This method gives me about 60-90 samples for each hour (60 samples if two years are available, or 90 samples if three years are available).

To give you an idea of the speed of computing on OSG, I have a job that would take about 2 hours to run on the meso4 server, which has 32 processors, and created 288 ~65MB files. In comparison, I submitted the same job to the OSG, which ran in about 15 minutes. The difference is that I ran the 288 jobs on meso4 in serial, whereas I ran the 288 jobs simultaneously on the OSG, each with between 8-48 processors, as soon as the computing resources become available. If you don't account for the time in the queue waiting for resources to become available, it takes less than one minute (typically between 10-60 seconds with an average time of 30 seconds) to complete the computations on the 8,000+ cores. The trade off is in the transfer time...transferring the files from the compute node (farmed out somewhere in the United States) to my OSG home directory, and then from the OSG home directory to meso4 via scp. The transfer of files from OSG to meso4 can take up to 24 minutes if done in serial (about 4-5 seconds per file, and need to transfer 288 files).



I am sifting through TONS of data. Each variable for each hour is about 1MB in GRIB2 format, that is downloaded from the Pando archive. Converting the field to NetCDF when we save the final results, bloats the file size. On the OSG I have had to request computers that have a minimum of 8GB of memory in order to complete the computations. In the future, when there are 2 or 3 years of additional HRRR data, this method of calculating statistics will require a computer with a minimum of 16GB of available memory. (Computing the max, min, and mean of the data set is not memory intensive. Computing percentiles is memory intensive, because all the values need to be stored and then sorted.)

Examples of data computed:


Note: the graphs shouldn't be interpreted as a continuum. Each month is independent from one another. This explains the steep steps in a value from month to month. These steps from month to month are caused because the maximums are the month maximum. For example, the maximum on the 23rd hour in May and the 0th hour in June are very different. This is because the maximum in May may be on on May 30th, while the maximum in June may be on June 30th. Thus, the line representing the maximum of each hour in each month cannot be considered a continuous.

Because I only have 2.5 years of data, a loop showing a map of the maximum or minimum values appears to show features of two underlying weather features move around the map. These are caused by the extreme events for the day


Most useful commands:
condor_submit name_of_job_file.submit
condor_q
condor_q -nobatch
condor_q blaylockbk
condor_rm 1234567890

The queue looks different whether you are in the normal command prompt, or if you have loaded a singularity image. I find that the queue within the singularity image is filled with someone else submitting jobs. Jobs I've submitted get through the queue really quick in the normal command prompt (when I'm not in the singularity image).
  • Data Transfer of large files is somewhat cumbersome. I use the scp command to move files form OSG to the CHPC environment.
  • I use the netCDF python module for writing the data, and use basic compression (comlevel=1) which turns a 180MB file to a 65MB file. In those files contains 9 statistical fields (max, min, mean, p1, p5, p10, p90, p95, p99) a field for latitude and longitude, and details on the number of cores used to compute the statistics, total time, and beginning and ending date.
  • I couldn't set up SSH keys on Putty, so I gave up.

Python multiprocessing is easy to utilize, but doesn't really gain me any ground for what I'm doing. Primarily, multiprocessing is used to speed up the downloads, but it doesn't appear it makes much difference in the compute time. Below shows a scatter plot for 288 files of the number of cores used in the computation in relation to the number of seconds it took to complete the computations. Something to keep in mind is that some of these are calculating statistics for a high sample size, while some are calculating for a fewer sample size.

Below is for the 2m Temperature file:

What do we learn about the number of samples in a computation? They can either run fast or slow for a given number of samples:



Running the same script on the OSG for another variable, this time 10 m wind, you get a different result. Here you see that more cores does increase the speed of the computations, to a point.

 Number of samples for each hour: April 18, 2015 is when the archive started.



I have plotted some of the data from my experiment:


Then, I can compare these statistics to the observed temperature and the HRRR analysis data:

Wind, with the 95th percentile



Thursday, May 18, 2017

Save a python variable using Numpy

This just might revolutionize the way I return to projects and test functions

You can save a python variable using numpy!!!

Sure wish I would have known this a long time ago. For example, if you have a large python dictionary, you can save the dictionary, and then load it and use it later.

Example:

SomeDictionary = {'Name': np.array([1,2,3,4,5]), 'DATES':np.array([1,2,3,4,5])}

np.save('PyDict.npy', SomeDictionary)

get_the_dict = np.load('PyDict.npy').item()

What a miracle :)

Thursday, April 13, 2017

Friday, March 24, 2017

Multiprocessing vs Multithreading

I need to download files fast. I want to swamp my network cables to download as much as it can consume. In the past I've used multiprocessing because it was easy to split jobs to several processors. But it is more efficient to use multithreading.

Multithreading uses one processor, but you can set up a queue so that it will continuously run as much as it can handle.

As an example, here is some code...

# Brian Blaylock
# March 24, 2017                           Yesterday Salt Lake had lots of rain

"""
Fast download of HRRR grib2 files with threading
"""

from queue import Queue 
from threading import Thread
from datetime import datetime, timedelta
import numpy as np
import urllib2
import re
import os

var='TMP:2 m'

def download(URL):
    # Code to download something based on a URL
    # ...
def worker():
    while True:
        item = q.get()
        print "number:", item
        download(item)
        q.task_done()

num_of_threads = 10

q = Queue()
for i in range(num_of_threads):
    t = Thread(target=worker)
    t.daemon = True
    t.start()

# List of URL's to download from
URL_list = ['...','...','...']

timer1 = datetime.now()
for item in URL_list:
    q.put(item)

q.join()       # block until all tasks are done


There is speed up when downloading, and the copper wires are saturated with as much as it can handle with about 10 threads or ten processors.
Here is a visual of what I learned:



Tuesday, February 14, 2017

not python, but I made a handy image viewer

This isn't Python, but I made a handy web page tool for viewing images in a directory. Just dump this photo_viewer.php page in a public_html web directory with images in it, and this page will let you view the images by hovering a mouse over buttons with the image names, clicking the buttons, or selecting images from an option box.

https://github.com/blaylockbk/Web-Homepage/blob/master/photo_viewer/photo_viewer.php

Example page: http://home.chpc.utah.edu/~u0553130/PhD/UWFPS_2017/time-height/photo_viewer.php


Thursday, January 5, 2017

Simple Basemap Example

You can do a lot with Python basemaps. Here is an example of a few things you can do...

Friday, November 11, 2016

Verifying GFS dewpoint data with MesoWest observations

For a class assignment, I needed to verify GFS model output with MesoWest data. I used pygrib to read in the grib file, then used some custom functions that grab data from the MesoWest API.

Here is the output figures and my code.


The output:







gist: