EGB13

Radar on Blue Marble image

This is a test run of the nation-wide (lower 48) weather radar animated on a Blue Marble image, with state outlines.  The YouTube link is here, and it’s embedded below, along with a short discussion of the necessary steps to make it, and what could make it better.

This is just about a day’s worth of animation from images archived at 10-minute intervals.  They are encoded at 29.97 frames per second to NTSC DVD standard with a 16:9 aspect ratio:

Not bad, though I don’t like those pulsing white splotches. I’ll discuss how those might be elliminated later. First come the steps to make what I have now.

Executive version

The earth does not change, at least as far as we’re concerned here. A section of the Blue Marble image corresoponding the area covered by the radar is cropped out of the original one time and saved as a GeoTIFF file. The state boundaries are static, too, so an image of those with the correct geographic extents is saved as a PNG image in which all but the lines are transparent. Each time a new radar image is downloaded, the earth image, the radar image (in a transparent GIF), and the state lines are layered into a single image in that order from the bottom up. At this point everything is still rectilinear, so the resulting composite image can be projected with a source SRS of EPSG:4326.

Detailed procedure

Radar image

NOAA makes a full-resolution composite of the individual national radar images available as a regularly updated GIF file from http://radar.weather.gov/Conus/RadarImg/ in the file latest_radaronly.gif, and its world file is latest_radaronly.gfw. The background of these images, where there is no radar reflection to show, is transparent in the GIF image. It merely needs to be layered on top of a suitable background. Its world file is:

0.017971305190311
0.000000000000000
0.000000000000000
-0.017971305190311
-127.620375523875420
50.406626367301044

The upper-left corner is at (rounded) 50.41 north latitude 127.62 west longitude. (Technically, that’s supposed to be the center of the pixel, but my arithmetic is not that fastitidous.) There are 0.0180 degrees per pixel, both of longitude and latitude. The size of the image is 3400×1600 pixels, so there are (3400 x 0.0180) = 61.20 degrees longitude and (1600 x 0.0180) = 28.80 degrees latitude in the image. (In practice I let the computer do full-precision math, but rounded numbers are much more presentable here.) The coordinates of the lower-right corner are at (127.62 – 61.20) = 66.42 E and (50.41 – 28.80) = 21.61 N.

Earth image

The NASA Earth Observatory’s Blue Marble project has a series of color images of the entire Earth’s surface composited from cloud-free satellite images.  I’ve used a 21600×10800-pixel image that has a resolution of 2km per pixel. That’s for the whole earth. I determined what section of that image corresponds exactly to the radar image by using the geometry values in the radar’s world file. That section is 3666×1725 pixels. I didn’t retain the crop parameters, so it’s left as an exercise to the reader if you want to do the same.  I used a command-line utility from ImageMagick to crop out that section — the image is way too big to load into a GUI program on any system that I own.

Plug into the world file calculator the geographic extent of the radar image as derived above along with the pixel size of the cropped US earth image, and out pops the world file for the US section of the Blue Marble image. Now that image can be used in GIS programs, including gdalwarp and Quantum GIS (QGIS).

State outline image

The radar image can be composited on the cropped US portion of the Blue Marble image, but it needs state outlines if you want to know more precisely where the weather is and where it’s headed.  Results showing proper alignment with those lines is also good for determining whether the raster images of the Earth and the radar are being properly projected.

QGIS will dynamically project vector layers to any projection it know about. Shapefiles of the US state boundaries are easy to come by — I got mine from the Census Tiger project. QGIS will also save a project as an image. Unfortunately, it does not save them as a GeoTIFF image; not my version, at least. In any case, I’m not using a GIS program to overlay the three images. So I need an image of the state outlines of exactly the same geographic extent and in the same spatial reference system as the radar and earth images. We need to make a transparent GIF or PNG of the state lines with the same rectangular “projection” and the same extents as the radar and earth images.

I opened the US earth image in QGIS with the project projection set to match the image (EPSG:4326) and dynamic projection enabled. QGIS does not dynamically project raster images (not my version, at least), so any raster images used in a QGIS project need to have the same projection, and the project projection must be set to match. Vector layers added to the project will then be dynamically projected in the same way and be placed in the proper position relative to the raster layers.

So a vector layer is added containing a shapefile of the US state boundaries. That is sized so that a bit of margin around the raster layer is showing and it is saved as a PNG image. The raster image is then turned off, leaving only the state boundaries showing, and that is saved again as a PNG.

Both images are the exact same size, extent, and location. They are opened in the GIMP as layers with the layer including the US earth image on top. The image is cropped to the exact size of the earth surface, the top layer with that earth surface is discarded, the white background in the remaining states-only layer is turned transparent, and the image is saved as a PNG with transparency. (The details of how to do all that is beyond the scope of this article.)

That transparent PNG image of the state outlines now has the exact same geographic extent as the US earth image and it has the same projection. Yay!

Putting it all together

We have three layers — earth, radar, states — but they’re all different sizes. Doesn’t matter. They’re the same extent and projection, and they have the same aspect ratio. They can simply be scaled to a uniform size and layered into a new image. Since we want to project that final image, though, we need to preserve or re-derive the geographic information. Because we are most interested in the radar image, and because we already have a world file for it, we’ll make the new image the same size as that. Using ImageMagick’s convert utility:

convert us-huge.jpg -resize 3400x1600 latest_radaronly.gif -composite statesoutline.png -composite -quality 95 usradar-huge-er.jpg

Hmm… that could be made more efficient.it was copied directly from the script I use from cron to do this all automatically. It resizes the US earth image (us-huge.jpg) to the size of the radar image every time. But it doesn’t scale the states outline image. I must have pre-scaled the states outlines, and it would make sense to do that for the US earth layer, saving the need to resize it over and over and over again as updated radar images are downloaded.

Regardless, the output of that command, usradar-huge-er.jpg, is an equirectangular (non)projection of all three layers combined. All that remains is to project it with gdalwarp. Its size, extent, and geometry match the radar image, so the world file that NOAA provides for that fits the new image, also. A copy of that in usradar-huge-er.jpgw is all that’s necessary for this to do the projection:

gdalwarp -s_srs EPSG:4326 usradar-huge-er.jpg -t_srs EPSG:2163 usradar-huge-laea.tif

EPSG id 2163 is the Lambert azimuthal equal area projection used by the National Atlas of the U.S. The output is a GeoTIFF file, though I’m not (currently) using this as input for any more GIS processing.

The last step is to add the timestamp. I get the string for that from the radar image file in the format I want (let me know if there’s an easier way for this):

tstamp=`perl -e 'my @f=stat("latest_radaronly.gif"); my @t=gmtime($f[9]); printf "%0.4d%0.2d%0.2d-%0.2d%0.2d",1900+$t[5],1+$t[4],$t[3],$t[2],$t[1];'`

And that is composited onto the projected image, writing the result to a new file that includes the timestamp in the filename:

convert usradar-huge-laea.tif -resize 1024x1024 -background black -fill yellow -font fixed -size 165x30 label:"${tstamp}" -geometry +$(((1024-165)/2))+6 -composite archive/usradar-${tstamp}.jpg

That includes resizing to a maximum 1024×1024 while preserving the aspect ratio, so the actual results in this case are 1024×597. The file goes into a subdirectory named “archive” where it is saved as a frame for the animation. Ten minutes later it happens again. Here is the whole script that cron invokes:

#!/bin/bash

URL="http://radar.weather.gov/Conus/RadarImg/latest_radaronly.gif"
DIR=/opt/mapdata/bluemarble

cd ${DIR} || exit 1

#- download
dlfile=latest_radaronly.gif
if [ -s ${dlfile} ] ; then
        since="-z ${dlfile}" #only download if it's newer than the last time
else
        since=""
fi
/usr/bin/curl $since --silent --remote-time -o ${dlfile} "${URL}"

#- conversions
convert us-huge.jpg -resize 3400x1600  latest_radaronly.gif -composite statesoutline.png -composite -quality 95 usradar-huge-er.jpg
gdalwarp -s_srs EPSG:4326 usradar-huge-er.jpg -t_srs EPSG:2163 usradar-huge-laea.tif
tstamp=`perl -e 'my @f=stat("latest_radaronly.gif"); my @t=gmtime($f[9]); printf "%0.4d%0.2d%0.2d-%0.2d%0.2d",1900+$t[5],1+$t[4],$t[3],$t[2],$t[1];'`
convert usradar-huge-laea.tif -resize 1024x1024 \
        -background black -fill yellow -font fixed -size 165x30 label:"${tstamp}" -geometry +$(((1024-165)/2))+6 -composite \
                archive/usradar-${tstamp}.jpg

Animation

I use ffmpeg to put the frames into an animation using this script:

#!/bin/bash
DIR=/opt/mapdata/bluemarble
cd ${DIR}/archive || exit 1
[ -d tmp.d ] || mkdir tmp.d || exit 1
n=1 ; ls | grep 'usradar-.*\.jpg' | while read f ; do
    [ -s $f ] && {
        ln $f tmp.d/t-${n}.jpg ;
        let n=$n+1 ;
    }
done
ffmpeg -r 29.97 -i tmp.d/t-%d.jpg -target ntsc-dvd -aspect 16:9 -an us-radar.mpg
rm -rf tmp.d

That creates links from the image files to sequentially numbered files in the tmp.d directory, and those are used as input to ffmpeg. It tests that the image files are not of zero size. That’s less reliable but faster than testing that they’re a valid image. That’s probably not an issue given how these images are created; it’s inherited from my other animation scripts.
The resulting animation has an aspect ratio of 16:9, which is 1.777778. The input frames are, by chance or design, 1024×597, which is 1.715243. That’s close enough and saves the expensive step of recompositing each frame onto a 16:9 matte.

Improvements

Those radar images, though only of the radar with the rest transparent, are obviously meant for a white background. That is why the areas of noise round the individual radar installations show up as white. In the downloaded transparent GIF, those areas are not transparent, and they’re almost, but not quite, white.
The cure that I see is expensive… for each pixel, convert the RGB colors to HSV, and then turn the pixel to transparent if the value is high and the saturation is low. “High” and “low” remain to be defined, probably empirically. That would have to be done prior to compositing the layers into the finished image. Right now I’m not ambitious enough to work out a program to do this.
I’m currently doing the compositing and warping each time the radar layer is downloaded and then overwriting the radar image the next time rather than preserving it. This method wouldn’t work with the composited images I’m archiving, so any animation I make from these I already have archived will be filled with pulsing white splotches.
The state lines are too thin when the image is scaled down for animation. I should make them thicker in the overlay file that contains them.

One last thought… GRASS GIS is a system that incorporates the GDAL library and a lot more. I haven’t wrapped my head around it yet, but it occurs to me that I should look into it in case it could be applied to some of this.