NOAA radar overlaid on a part of NASA's Blue Marble image of the southeast US.
The gdal set of GIS utilities is nothing new, and they are even incorporated into other GIS programs. But that does not mean that everyone who will ever need to master them has already done so, and it does not mean that mastering them is simple. It’s not really difficult to start using utlities like gdalwarp directly, but there are details that keep it from being really simple. This article is part of my notes to myself as I try to learn some of this.
I wrote an article about using gdalwarp to perform a pair of standard projections from a rectilinear view of the Earth using a Blue Marble image from NASA’s Earth Observatory. Below I discuss using a custom projection to do something marginally useful — combining the Blue Marble with some NOAA weather imagery.
Below are the steps that went into making that image. Note that I work in Linux, but Windows and likely Mac versions of these tools exist.
The radar image
A transparent GIF radar-only image from NOAA
NOAA has a large scale mosaic image of weather radar for the 48 contiguous states available. It is also available in sections, and and those sections also are available with the radar portion only. Those have no state lines, highways, land or ocean, and the portion not covered by the radar image is transparent The images are GIFs with transparency. NOAA also makes a world file available for the individual sections, and we know how useful those are when using gdalwarp.
The southeast had the most interesting weather at the time I started working on this, so that is the section I used. The radar-only sector image is named southeast_radaronly.gif, and the world file for it is southeast_radaronly.gfw, which looks like this:
0.017971305190311 0.000000000000000 0.000000000000000 -0.017971305190311 -90.24006072802854 36.928147474567794
That tells us that the upper left corner is at roughly 90.24° west 36.93° north, and each pixel is 0.0180° in both width and height.
The pixel dimensions of the image are 840x by 800y, so the geographic width is 840×0.0180°=15.0960° and the height is 800×0.0180°=14.377°. Therefore the geographic coordinates of the lower left corner (not rounding) is determined by:
(90.24006072802854°-15.095896°) W (36.928147474567794°-14.377044°) N = 75.144165° W 22.551103° N
The Blue Marble land image
A larger (not the largest!) Blue Marble rectilinear image of the whole earth is at http://veimages.gsfc.nasa.gov/2433/land_shallow_topo_21600.tif. I didn’t make that a clickable link because it’s likely way too big for your browser to chew on. I downloaded it and worked on it with command-line utilities. Its pixel dimensions are 21600 by 10800. The upper left corner is at 180 W 90 N, and it is 360° of longitude wide and 180° of latitude tall. Using the world file calculator we generate this world file for the big Blue Marble image:
0.016666666666666666 0.00000 0.00000 -0.016666666666666666 -180.0 90.0
So we know that each pixel is 0.01667° (rounded) tall and wide.
Next we determine the pixel-location within the Blue Marble image of location 90.24006072802854° W 36.928147474567794° N, which corresponds to the upper left corner of the radar image:
(180-90.24006072802854)/0.01666666666666667 = 5385.596356 -> 5386 x
(90-36.928147474567794)/0.01666666666666667 = 3184.311152 -> 3184 y
What is the Blue Marble pixel width & height of the area corresponding to the southeast radar tile?
(15.095896 wide)/0.01666666666666667 = 905.748000 -> 906 px wide°
(14.377044 high)/0.01666666666666667 = 862.622640 -> 823 px high
A portion of NASA's Blue Marble image showing the southeast US
That is all that is needed to crop a tile from the Blue Marble image that corresponds to the area of the radar image. We can use convert from the ImageMagick toolbox to extract that tile:
convert world-huge.jpg -crop 906x863+5386+3184 southeast-marble.jpg
Those same geographic coordinates and pixel sizes go into the world file calculator to create this world file for the tile cropped from the Blue Marble image:
0.01666666666666667
0.00000
0.00000
-0.01666666666666667
-90.24006072802854
36.928147474567794
Since we created the tile image with the name southeast-marble.jpg, these values go into a file called southeast-marble.jpgw (or .jgw).
Note that the scale of the pixels is the same as for the Blue Marble image, since we didn’t resize anything, and the geographic coordinates are identical to those of the radar image that we got from NOAA. We could have done that without the calculator.
Image projection
In order to warp both the radar image and the tile cropped from the Blue Marble into an Albers equal area projection centered on that area, we choose the central meridian to be the longitude of the center of the tile, and the two standard parallels to be located at 1/3 and 2/3 of the way through the image vertically:
Longitude = 75.144165°+(15.095896°/2) = 82.692113° W Parallel 1 = 22.551103°+(14.377044°/3) = 27.343451° N Parallel 2 = 22.551103°+(14.377044°/3)*2 = 32.135799° N
The proj4 notation for that projection is:
+proj=aea +lon_0=82.692113w +lat_1=27.343451n +lat_2=32.135799n +datum=NAD27 +ellps=clrk66 +units=m +no_defs
Note: The datum, ellps, units, and no_defs parameters are not mentioned for the aea projection in the manual for proj version 3. (I didn’t look at the version 4 manual.) I took those from an aea projection definition in the Quantum GIS program config dialog. (Gdalwarp works without them, but QGIS, used in a later step, required them, so I used the same notation in each program.) I could have chosen a different datum and ellipse, but I don’t know that it would matter too much at these scales. It certainly doesn’t matter to me for the purpose I have in mind.
Now we can project the radar and background tiles:
gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aea +lon_0=82.692113w +lat_1=27.343451n +lat_2=32.135799n +datum=NAD27 +ellps=clrk66 +units=m +no_defs" southeast_radaronly.gif radar.tif
gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aea +lon_0=82.692113w +lat_1=27.343451n +lat_2=32.135799n +datum=NAD27 +ellps=clrk66 +units=m +no_defs" southeast-marble.jpg land.tif
Projected radar-only
Southest Blue Marble image in AEA projection
Put it all together
That gives us two images of the same area, one radar and one of the land surface, both warped with a common custom projection centered on the area. They were both created as GeoTIFF files by gdalwarp, so their projection information is included in the image file and separate world files for these are not needed.
At this point I imported them into QGIS. I’m not going to include a tutorial of that here, but an outline of what I did is this…
I created a custom projection in QGIS identical to the target SRS above. I set the project default projection to that and enabled on-the-fly projection.
I imported both the land and the radar projections as raster layers into QGIS and stumbled upon how to enable the transparency in the radar image. They seemed to correspond with each other, which increased my confidence that everthing was done correctly thus far.
Next I added the state outlines from a Tiger shapefile. The coastlines aligned with the land surface layer, indicating that the Blue Marble tile was cropped and projected properly. And the locations of the storms relative to the state lines jibes with the composite image from NOAA. Everything is where it is supposed to be. Success!
NOAA radar overlaid on a part of NASA's Blue Marble image of the southeast US.
Except for the QGIS part this could all be automated in a script. An alternative procedure might be to scale the Blue Marble tile to the same size as the radar sector before projection, overlay the radar on the land image strictly as raster images of equal dimension, and then to warp them. That still leaves out the superposition of the state lines. I haven’t got around to working on that yet.
