Meteorological meanderings, part 2: visualizations and long-range forecasts

In the previous post, we used GDAL to query temperature information from a downloaded HRRR GRIB2 file. That query, however, was for a single point. Sometimes it helps to be able to see many points at once, preferably on a map. This can be useful for determining whether a particular band “looks right”, as well as for doing sanity checks when writing code to query or visualize the data yourself.

A brief look at QGIS

The best tool I’ve found so far for visualizing meteorological data is QGIS. It’s an extremely powerful cross-platform piece of software, with many features, but we don’t need the vast majority of them. To visualize a band from HRRR with QGIS, we can start by bringing in a world map to orient ourselves. Thankfully, that’s just a few clicks away, thanks to built-in support for OpenStreetMap:

Animated screen capture of loading OpenStreetMap in QGIS

Once we have the map pulled up, we can load HRRR as a separate layer on top of the OSM layer:

Animated screen capture of loading HRRR data in QGIS

The HRRR file appears twice in the file browser. That is because it can be loaded as either a normal raster file, or as a mesh file. Loading it as a raster file is much faster, but it’s more difficult to control the appearance of layers, especially temporal data, which is why the second file — the mesh — is chosen above. After the data is loaded, which can take a while (that part isn’t shown), we can go through its properties to select the band we want to see:

Animated screen capture of picking the Temperature band in QGIS

Setting the opacity to 50%, as seen above, helps to retain that orientation given to us by the OSM layer. After we have the desired band showing, we can use the “Identify Features” mode in QGIS to get details about the band at specific locations:

Animated screen capture of accessing a single data point from the Temperature band in QGIS

In this particular data file, it’s approximately -4 degrees Celsius in the vicinity of Indianapolis, Indiana. This is a very cursory look at QGIS, but it’s enough to show how useful it can be for visualizing meteorological data.

Long-range forecasts and their projections

To get a long-range forecast, GFS is a useful model. It’s global, and it goes out to more days in the future than is reasonable. GFS data, like that of HRRR, is also downloadable from NOMADS. Accessing it via GDAL, however, poses a slight problem. Armed with knowledge that temperature at 2m above ground level is in band 435, we should be able to do this query just like we did with HRRR:

$ gdallocationinfo -b 435 -wgs84 gfs.t00z.pgrb2.0p25.f192 -86.1583502 39.7683331

Report:
  Location: (-345P,201L)

Location is off this file! No further details to report.

Oops! That’s not good. We’re using the same coordinates, specified as WGS84. GFS is supposed to have data for the entire globe, so how can the specified location be off the grid? Before digging into that, let’s extract the temperature band into a separate file, so that there’s less data to deal with. To do that, we can use another GDAL tool, gdal_translate. Its primary purpose is translating between different raster data formats, but it can also be used to extract bands. This is how we can use it to extract a single band and, at the same time, convert it to GeoTIFF:

$ gdal_translate -b 435 gfs.t00z.pgrb2.0p25.f192 gfs.t00z.pgrb2.0p25.f192-temperature.tif

Input file size is 1440, 721
0...10...20...30...40...50...60...70...80...90...100 - done.

Now let’s see what gdalinfo can tell us about this band from its own GeoTIFF file:

Driver: GTiff/GeoTIFF
Files: gfs.t00z.pgrb2.0p25.f192-temperature.tif
Size is 1440, 721
Coordinate System is:
GEOGCRS["Coordinate System imported from GRIB file",
    DATUM["unnamed",
        ELLIPSOID["Sphere",6371229,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (-0.125000000000000,90.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -0.1250000,  90.1250000) (  0d 7'30.00"W, 90d 7'30.00"N)
Lower Left  (  -0.1250000, -90.1250000) (  0d 7'30.00"W, 90d 7'30.00"S)
Upper Right (     359.875,      90.125) (359d52'30.00"E, 90d 7'30.00"N)
Lower Right (     359.875,     -90.125) (359d52'30.00"E, 90d 7'30.00"S)
Center      ( 179.8750000,   0.0000000) (179d52'30.00"E,  0d 0' 0.01"N)
Band 1 Block=1440x1 Type=Float64, ColorInterp=Gray
  Description = 2[m] HTGL="Specified height level above ground"
  Metadata:
    GRIB_COMMENT=Temperature [C]
    GRIB_DISCIPLINE=0(Meteorological)
    GRIB_ELEMENT=TMP
    GRIB_FORECAST_SECONDS=691200 sec
    GRIB_IDS=CENTER=7(US-NCEP) SUBCENTER=0 MASTER_TABLE=2 LOCAL_TABLE=1 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2021-01-15T00:00:00Z PROD_STATUS=0(Operational) TYPE=1(Forecast)
    GRIB_PDS_PDTN=0
    GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES=0 0 2 0 96 0 0 1 192 103 0 2 255 0 0
    GRIB_PDS_TEMPLATE_NUMBERS=0 0 2 0 96 0 0 0 1 0 0 0 192 103 0 0 0 0 2 255 0 0 0 0 0
    GRIB_REF_TIME=1610668800 sec UTC
    GRIB_SHORT_NAME=2-HTGL
    GRIB_UNIT=[C]
    GRIB_VALID_TIME=1611360000 sec UTC

Yep, it’s certainly the band we wanted. But above the band-specific information, there is a clue as to what’s wrong. The reported corner coordinates are, approximately, (0, 90) at the upper left and (360, -90) at the lower right. While the Y-coordinates seem fine for latitude, the X-coordinates are not. We want to see longitude represented as -180 to 180, not 0 to 360. That explains why the location information request failed: -86.1583502 really is off the grid here!

To fix this issue, we can use GDAL’s warping capability:

$ gdalwarp -t_srs "EPSG:4326" -wo SOURCE_EXTRA=100 --config CENTER_LONG 0 gfs.t00z.pgrb2.0p25.f192-temperature.tif gfs.t00z.pgrb2.0p25.f192-temperature-warped.tif

Creating output file that is 1440P x 721L.
Processing gfs.t00z.pgrb2.0p25.f192-temperature.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

The above command does a few things. First, it specifies a different coordinate reference system, or spatial reference system (SRS). EPSG:4326 is another name for WGS84. Next, it grabs a bit of extra data for the warp procedure, to ensure no gaps exist when stitching. Finally, the config option CENTER_LONG 0 is what tells GDAL to recenter the image, which will get us the coordinate range we want:

$ gdalinfo gfs.t00z.pgrb2.0p25.f192-temperature-warped.tif

Driver: GTiff/GeoTIFF
Files: gfs.t00z.pgrb2.0p25.f192-temperature-warped.tif
Size is 1440, 721
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-179.999755859375000,90.125000000000000)
Pixel Size = (0.249999673885899,-0.249999847496904)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-179.9997559,  90.1250000) (179d59'59.12"W, 90d 7'30.00"N)
Lower Left  (-179.9997559, -90.1248900) (179d59'59.12"W, 90d 7'29.60"S)
Upper Right ( 179.9997745,  90.1250000) (179d59'59.19"E, 90d 7'30.00"N)
Lower Right ( 179.9997745, -90.1248900) (179d59'59.19"E, 90d 7'29.60"S)
Center      (   0.0000093,   0.0000550) (  0d 0' 0.03"E,  0d 0' 0.20"N)
Band 1 Block=1440x1 Type=Float64, ColorInterp=Gray
  Description = 2[m] HTGL="Specified height level above ground"
  Metadata:
    GRIB_COMMENT=Temperature [C]
    GRIB_DISCIPLINE=0(Meteorological)
    GRIB_ELEMENT=TMP
    GRIB_FORECAST_SECONDS=691200 sec
    GRIB_IDS=CENTER=7(US-NCEP) SUBCENTER=0 MASTER_TABLE=2 LOCAL_TABLE=1 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2021-01-15T00:00:00Z PROD_STATUS=0(Operational) TYPE=1(Forecast)
    GRIB_PDS_PDTN=0
    GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES=0 0 2 0 96 0 0 1 192 103 0 2 255 0 0
    GRIB_PDS_TEMPLATE_NUMBERS=0 0 2 0 96 0 0 0 1 0 0 0 192 103 0 0 0 0 2 255 0 0 0 0 0
    GRIB_REF_TIME=1610668800 sec UTC
    GRIB_SHORT_NAME=2-HTGL
    GRIB_UNIT=[C]
    GRIB_VALID_TIME=1611360000 sec UTC

That looks much better. The upper left coordinate is now (-180, 90) and lower right is (180, -90). This is what we expect to see. Now that we’ve confirmed that the corner coordinates look good, we can perform our temperature query against the warped file:

$ gdallocationinfo -wgs84 gfs.t00z.pgrb2.0p25.f192-temperature-warped.tif -86.1583502 39.7683331

Report:
  Location: (375P,201L)
  Band 1:
    Value: -1.61920776367185

Hooray, it works! Using this procedure, we can convert any band to an easily queryable format. And, just like with HRRR, if we bring in more bands, we should be able to generate a multi-day forecast from GFS data.

We’re not going to cover this now, but there’s another interesting conversion process we could get involved in, if we wanted to access global forecasting data from Deutscher Wetterdienst, the German Meteorological Service. They produce the ICON model, which uses a triangle-based grid system. Thankfully, they have an extraordinarily well-written document describing how to convert that data to the familiar GRIB2 format.

Next up, we’ll switch gears and look at performing geocoding and reverse geocoding using a locally-hosted service.