Category Archives: NCL

Displaying WRF data using OpenLayers

Introduction

So, you’re bored of well made, clear, beautiful but static plots which can show people in one glance what the weather is doing. You want exciting slippy maps that people can zoom around and go “oooooohh pretty!”. I’ve attempted this several times in the past, and each time been thoroughly beaten, as my lovely pressure field from WRF ends up as a single pixel somewhere off the coast of east Africa.

Well the good news is, there is a comparatively easy method of doing this, thanks to a bunch of libraries and programs which do the hard work. The even better news is that if you are already producing plots using NCL, then it just requires a couple of modifications to your scripts.

My downfall in the past was always re-projecting WRF data. Because WRF is a grid based model, where the curved atmosphere is represented by rectangular arrays, your WRF data is already projected. You have hopefully chosen your projection to minimise distortion over your domain, and unless you have very specifically designed your WRF domain for web mapping, this won’t be the same system used by common web tools, and you will have to re-project your data before you can overlay it on a web map. There are a number of points where you could do the re-projection:

  1. The raw netcdf files can be re-projected using a post-processing tool such as UPP
    • UPP is a bit fiddly, and returns grib files
  2. The arrays within the netcdf files can be re-projected using GDAL,NCL, or Python
    • handling the variables within the netCDF files is tricky with GDAL ( see Reading netcdf files with GDAL and re-project files with GDAL for hints)
    • doing re-projection well in Python is a bit of a pain
    • NCL is probably the most adept at this, particularly with the Earth System Model regridding tools in latest versions
  3. The images produced by NCL or other tools can be made on a new projection
    • NCL or Basemap in Python can re-project at the time a plot is produced
    • The image needs to take up the whole frame, without extra stuff
  4. The images produced by NCL can be re-projected after the event using GDAL

I find Option 4 by far the easiest way, since the excellent VAPOR package lets you export properly georeferenced geotiffs directly from NCL.

Dependencies

  • NCL
  • VAPOR
  • GDAL
  • ImageMagik

Geotiffs from NCL and VAPOR

Vapor includes the script $VAPOR_ROOT/share/examples/NCL/wrf2geotiff.ncl. See the examples here for full instructions.

I had to set the workstation type to oldps, not ps as stated in the docs, and I also had to output one frame per geotiff in order for the re-projection to work. The key modifications to an NCL script are:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

; load the wrf2geotiff library from the VAPOR distribution:
load "$VAPOR_ROOT/share/examples/NCL/wrf2geotiff.ncl"

; do everything else as normal, e.g. get data, loop over time, etc. etc.
....

type      = "oldps"
wks       = gsn_open_wks(type,plot_name) ; I open and delete a workstation each timestep
wrf2gtiff = wrf2geotiff_open(wks)

; set your plot, map, contour options etc
mpres = True

; make the plots
contour_psl = wrf_contour(a,wks,slp,opts_psl)

; overlay on the map
plot = wrf_map_overlays(a,wks,(/contour_psl/),pltres,mpres)

; write into the geotiff object
wrf2geotiff_write(wrf2gtiff, a, times(it), wks, plot, True)

; advance frame between time steps after the wrf2geotiff_write
frame(wks)

; this triggers actual writing of all geotiff output frames
wrf2geotiff_close(wrf2gtiff,wks)

delete(wks)
delete(wrf2gtiff)

Reprojecting: the easy way

Assume we now have the geotiff precip.tiff. Those clever folks at vapor have correctly georeferenced it and set all the metadata. I’m using a Lambert Conformal Conic projection, and I can look at the projection info using gdalinfo.

$>gdalinfo precip.tiff

Size is 995, 659
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT[,0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",59.5],
    PARAMETER["standard_parallel_2",49.5],
    PARAMETER["latitude_of_origin",56],
    PARAMETER["central_meridian",-1.5],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["unknown",1]]
Origin = (-1998307.075827066786587,1323201.405077444389462)
Pixel Size = (4020.745950011026252,-4021.886754136362924)
...

This is very useful info, even if you want to try a different approach, since sometimes it’s hard to work out how the WRF projection options translate into WKT or Proj4 strings.

The only concerning thing is that the spheroid is not set, and defaults to WGS spheroid, whereas WRF uses a sphere. Perhaps this might be fixed in a later version, or you could fix it with GDAL to set the datum. It could lead to features being displaced, see this paper. However, for coarse display purposes I’ll ignore that for now.

The “Web Mercator” projection

The most commonly used projection in the web mapping world (google maps, bing, osm) has become known as the web mercartor. There is an excellent article about this projection. The short version is that it now has the code EPSG:3857. The downside is that many conversion routines don’t handle it properly, since it specifies points on the earth using the WGS1984 spheroid, but projects these using a spherical earth.

It easy to re-project the geotiff to web mercartor projection using gdalwarp:

$> gdalwarp -t_srs EPSG:3857 precip.tiff  precip_warped.tiff

Check it looks right, then also grab some info which we will need later:

$> gdalinfo precip_warped.tiff
Size is 1322, 780
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
AUTHORITY["EPSG","3857"]]
Origin = (-4708884.528225043788552,10384546.176803732290864)
Pixel Size = (6876.374362466389357,-6876.374362466389357)

Corner Coordinates:
Upper Left  (-1998307.076, 1323201.405) ( 42d18' 2.27"W, 62d31' 7.21"N)
Lower Left  (-1998307.076,-1327221.966) ( 25d12'38.84"W, 41d 3'55.88"N)
Upper Right ( 2002335.144, 1323201.405) ( 39d21'56.03"E, 62d29'56.26"N)
Lower Right ( 2002335.144,-1327221.966) ( 22d15'18.11"E, 41d 3'13.56"N)
Center      (    2014.034,   -2010.280) (  1d28' 3.44"W, 55d58'54.76"N)

Now I recode the no-data values as something which will become transparent. This may depend on what you are plotting, but for precipitation, we can safely convert no-data to 0.

$>gdal_translate -a_nodata 0 precip_warped.tiff precip_mask.tiff 

Then I convert to a png, with what was encoded as nodata as set to transparent thanks

$> gdal_translate -of PNG -scale precip_mask.tiff precip.png   

Finally, I needed to add a hack. Although I tried to make my lowest contours (i.e. no rainfall) transparent, they came out white after all that conversion. So I also convert white to transparent using ImageMagik, so the rain looks better overlaid on a map:

$>convert precip.png -transparent white precip_final.png

I wrap up those commands in a bash script:

#!/bin/bash

for f in "$@"
do
    bname=`basename "$f" .tiff`
    echo $bname
    gdalwarp -t_srs EPSG:3857 $f $bname.tmp1.tiff
    gdal_translate -a_nodata 0 $bname.tmp1.tiff $bname.tmp2.tiff
    gdal_translate -of PNG -scale $bname.tmp2.tiff $bname.tmp3.png
    convert $bname.tmp3.png -transparent white $bname.png
    rm *.tmp?.tiff
    rm *.tmp?.png
done

Overlay in OpenLayers

Once you have the re-projected png, it is pretty straightforward to use it as an image overlay. The size and bounds can be got from the output of gdalinfo as above. I’ll do something fancier in a later post, once I have some nice animation controls, but the important javascript snippets look something like:

<script src="../OpenLayers/OpenLayers.js"></script>

  function init() {
    map = new OpenLayers.Map("mapCanvas");
    var osm            = new OpenLayers.Layer.OSM();

    var imgOpt = {
      isBaseLayer: false,
      maxResolution: "auto",
      resolutions: map.layers[0].resolutions,
      projection: map.getProjectionObject(),
      strategies: [new OpenLayers.Strategy.Fixed()],
      displayInLayerSwitcher: true,
      transparent:true, 
      format:"image/png", 
      alpha:true, 
      opacity:0.4 
    };
    var imgBounds =  new OpenLayers.Bounds(-4708884.528, 5020974.174, 4381682.379,10384546.177);
    var imgSize   = new OpenLayers.Size(1322, 780);
    imgUrl        = "precip_final.png";        
    image         = new OpenLayers.Layer.Image("precip", imgUrl,imgBounds, imgSize,imgOpt);

};

WRF in Google Earth

Traditionally weather maps aimed to present an enormous array of data as a simple output: a sun, a sun behind a cloud, or some rain.  Or possibly a cloud getting blown along by a face.  That’s great, as what most people want to know is: will it rain or not?

But this hides is the beauty and complexity behind the weather, we are immersed in a big mass of swirling fluid, rising, falling, heating, cooling.  Enormous amounts of energy are transported around the world, air masses collide, jet streams arc and split. In the midst of this, occasionally everything stills, and there is hardly a breeze.  It’s one big never-ending dance, and it looks beautiful.

That’s why I’m trying to put forecast data into Google Earth, because you get a real feel for the fact the UK’s weather doesn’t end at the edge of the map.

Using NCL and some code from rpavlick, this is my first attempt to put a forecast into Google Earth. At the moment it is just a picture, but hopefully I’ll get a Google Earth gadget working soon.