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

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' % service

image = cimgt.GoogleTiles(url=url)

ax.add_image(tile, 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