Month: May 2019

Originally published 2009-03-26.

Long story - short

To project the individual frames for this video I used:

gdalwarp -tr 10000 10000 -s_srs EPSG:4326 -t_srs EPSG:2163 nation.gif nationproj.tif

Long story - long

I finally figured out some of the mojo necessary to get some good out of gdalwarp. For me the key was to understand the "spatial reference system" (SRS) parameters. Key to understanding that was grokking the "+proj" notation used in the proj program.

It helped that I already had a handle on some of the basics of map projections. Geek that I am, one of the titles on my light reading list is Flattening the Earth: Two Thousand Years of Map Projections by John P. Snyder, who appears to be the go-to man for map projections. Along with interesting discussion of the pros and cons of each projection, he also gives the formulas for the forward projection of each. Using those formulas (and perlMagick) I wrote my own working but painfully slow perl programs to implement a couple of conic projections. But with almost fifty thousand images to project those programs proved entirely impractical, if educational. I had plenty of time to dive into gdalwarp while my programs churned away for weeks on end.

That trivia aside, here is what you need to know, or what I needed to know, to get started with gdalwarp. Warning: I am not a professional geographer or cartographer or anything of the sort. My use and understanding of some of the technical terms may make such a professional wince. But we don't care because we just want to make the software work for us, regardless of how fractured the information bouncing around inside our heads is.

Spatial Reference systems & projection parameters

Gdalwarp uses a source SRS to interpret the geometry of the source image, and it uses a target SRS to know about how you want to reproject it. These can be specified in "proj" format as described in the documentation for the proj program. At first glance that stuff looks a little like heiroglyphics, but a little knowledge of map projections helps to bring it our of the fog. For example, the formulas for the Albers Equal Area projection require as constants a central meridian and two standard parallels of latitude. Both the proj program and gdalwarp have the Albers Equal Area projection, among many others, built in. You can specify it with "+proj=aea" followed by the parameters it needs. E.g., "+proj=aea +lon_0=90w +lat_1=20n +lat_2=60n" for a projection centered on 90 degrees west longitude with standard parallels 20 and 60 north.

Fortunately, a shorthand notation is available for common projections. You do not have to remember and parrot the parameters for the Lambert Azimuthal Equal Area projection used by the National Atlas of the U.S. ("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), or know that the Cassini projection used for the Kertau 1968 Singapore Grid uses "+proj=cass +lat_0=1.287646666666667 +lon_0=103.8530022222222 +x_0=30000 +y_0=30000 +a=6377304.063 +b=6356103.038993155 +towgs84=-11,851,5,0,0,0,0 +units=m +no_defs". Those parameters are built into gdalwarp under the handy EPSG ids of 2163 and 24500, respectively. (I assume that second one is, though I haven't tried it… I'd certainly like it available simply as EPSG:24500 if I needed it.)

You do need a source for those EPSG ids. I find the projection setup in Quantum GIS to be one handy source since it's already on my desktop. That's where I found the Singapore example. Online sources include the EPSG Geodetic Parameter Registry and SpatialReference.org

World files

There is one more bit of information about the source image that gdalwarp needs to do its thing. There must be some way for the program to determine for each source image pixel location the corresponing x and y value in the source SRS -- that the pixel at [112,14] is for 32.25N 105.10W, e.g. If you have a GeoTIFF file as your source, you're in luck, because all the information needed is there. If it's a GIF or PNG or JPG file, or a non-Geo TIFF file, there is fortunately a way to provide that information, if you have it, in the form of a world file. A world file has the same base name as the image file and a .suffix of .gfw, jfw, .pgw or tfw for a GIF, JPG, PNG or TIFF file, respectively.

World files won't work for input images in which the magnitude of the x and y values differ across the image. It fits the bill for rectangular latitude/longitude non-projected images, and it's perfect for UTM. The format of a world file is simple and is described in that Wikipedia article. The situation is rather trivial if your input file is a simple non-rotated grid of equally spaced latitude and longitude coordinates, like this image:

That is an unprojected image with a constant change of longitude for each x pixel and a constant change of lattitude for each y pixel. If you know the coordinates of the upper-left corner of the image, if it's not rotated, and if you know its pixel dimensions, and if you remember the arithmetic you learned in grade school, you can easily figure out the values for the world file. (If it's rotated, study that Wikipedia world file article for me.)

When I was using my homemade projection program I guestimated the latitude and longitude of opposing corners of that radar image. I didn't realize then that the NWS has given us exact values in the form of...a world file! In this directory listing on the NWS website there are world files for many of the images they supply there. They do not supply the value for the small image that I'm using, but the lat-lon of the corners is the same as for the large image, so simple arithmetic gives me the degrees per pixel values I need to change in the world file. It comes out looking like this:

0.10183700000000
0.000000000000000
0.000000000000000
-0.101965000000000
-127.620375523875420
50.406626367301044

Update: I have created an online calculator for world files.

Putting it all together

In order to use the simple non-projected lat-lon format of these images you need to specify a corresponding source SRS for gdalwarp, and fortunately one is defined as EPSG id 4326. [Update: see this link, which I'm still trying to digest.] So we're getting close to finishing that gdalwarp command line at the top of this page. Here is the explanation of each part of the command:

Parameter	Significance
-tr 10000 10000	Hmm. I knew you were going to ask about this one. The gdalwarp manpage describes this as "output file resolution (in target georeferenced units)". It determines the size of the output, and I found a reasonable set of values to use purely by trial and error.  Further experimentation shows that you can leave it off and gdalwarp will use reasonable defaults.
-s_srs EPSG:4326	Source is a lat-lon grid of pixels
-t_srs EPSG:2163	The destination image has the same projection as the National Atlas for the resulting image. If it's good enough for them"
nation.gif	My input image, from the NWS website, is a copy of latest_Small.gif.
nationproj.tif	The target file, which will be created as one of those handy GeoTIFF files that other GIS programs can use.

That's all it took for me to warp that national radar image into a reasonable projection. In hindsight I should have used an equidistant projection instead of the equal area, so that the speed of the storms might be more constant across the whole continent. But at this scale it probably doesn't matter too much:

I converted the TIFF to PNG for the web, but gdalwarp can write a GIF file directly if you want. By the way, you can tell that this image is wider than what you see in the video. This is the correct projection. When I created the original MPEG file I preserved this shape in a 16:9 aspect ratio frame, same as HDTV, but it looks like after all this YouTube mashed it into a 4:3 area. That's a topic of investigation for another day.

I have used gdalwarp on other images. I've edited a GeoTIFF image in the GIMP, which does not preserve the georeferencing tags. I've added them back with gdalwarp: The listgeo program can create a world file from the original GeoTIFF, and that world file can then be used to reproject the edited file back into a GeoTIFF (using the same SRS for source and destination). Both listgeo and gdalinfo, by the way, can give you information about a GeoTIFF file if you need it, including the EPSG ids of the projection in use.

I'll finish this overlong treatise with a warning: gdalwarp is picky about the format of the TIFF files it will read. Any GeoTIFFs you find ought to already be of the proper form. Otherwise it sometimes likes to complain in cryptic ways about bands and such. In such a case, Google is your friend, even if gdalwarp refuses to be.

Final note

I'm not an expert! If you made it this far you now know practically as much as I do on the subject of gdalwarp. If this leaves you with further questions, remember that Google really is your friend.