A Gentle Introduction to GDAL, Part 2: Map Projections & gdalwarp
(If you’re new to GDAL, you might want to start with Part 1, which covers installation,
From ground level, an unobscured view of the horizon looks ruler-flat, an illusion that persists even from the roughly 30,000-foot-altitude of an airliner or the summit of Mount Everest. It’s only at twice that height, an elevation few of us will ever reach, that the first hint of a curve becomes visible.
This view was first seen in 1935 by two Army officers floating above the plains of South Dakota, aboard the Explorer II helium balloon. From a press release issued by National Geographic:
From its vantage point 72,395 feet in the air, the highest point ever reached by man, the camera registers the horizon 330 miles away, sweeping like a great arch across the photograph. The straight black line ruled across the top brings out the curvature of the Earth.
A balloon can fly into the stratosphere, but no higher, so it wasn’t until the advent of the space race that we got a better view of the edge of the Earth. From an altitude of several hundred miles one’s view stretches for thousands of miles, and the curvature of the Earth is clear.
From an even higher vantage point (900,000 miles), the Earth becomes a flat disk on a black field. (Even at this distance the stars are likely washed out by glare from the Earth.)
Each of these views is a distorted, two-dimensional representation of a three-dimensional surface. The appearance of the Earth is determined by the vantage point and optical characteristics of the camera taking the picture. Some features (those directly under the observer) are preserved very well, but others become quite deformed—it’s impossible to take a photograph that accurately represents the distances, areas, and angles of every part of a 3D object.
This is the fundamental challenge faced by cartographers—how to represent the 3D surface of the Earth on a flat piece of paper or computer screen? It’s also a challenge that mapmakers faced long before anyone ever saw the horizon as a curve.
Since a perfect representation isn’t possible, cartographers have developed an incredible variety of map projections, each designed to solve a particular problem. Every map is a compromise favoring one or more projection properties: area, form (or angle), distance, and direction.
Do you need to navigate a ship across the Atlantic? Use the Mercator projection, which keeps lines of latitude parallel. Show global temperature? Use a projection that preserves area, like Mollweide. Make a reference map? Use a compromise projection like Robinson. Show Antarctica or sea ice? Use polar stereographic.
To get started, download a global geo-referenced image like Natural Earth I with Shaded Relief and Water (zip file). This is a very nice map of the world with muted colors, based on NASA’s Blue Marble (yeah!) with shading for mountain ranges and the ocean floor.
Once you’ve downloaded and unzipped the TIFF, run
gdalinfo NE1_50M_SR_W.tif to take a look at the metadata, particularly the coordinate system:
Coordinate System is: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"]]
Of particular note are the units
degrees and the line starting with “authority”
"EPSG","4326". Unit indicates that each pixel is related to a real-world unit, not just pixels—degrees of latitude and longitude. The bracketed entries after authority (detailed explanation) are a shorthand reference for the map projection—EPSG:4326—which is simply a grid of equal latitude and longitude. This results in a world map that has a 2:1 aspect ratio, and goes from 180˚ west to 180˚ east and 90˚ north to 90˚ south. Of course it’s a bit more complicated than that (in ways I’ll describe later) and the full definition is lengthy, so having a simple 4-digit number makes things easier to type.
The file is big (10,800- by 5,400-pixels), so let’s make it a bit smaller using gdal_translate:
gdal_translate -r lanczos -tr 0.1 0.1 -co COMPRESS=LZW NE1_50M_SR_W.tif NE1_50M_SR_W_tenth.tif
If you read Part 1, most of this should look familiar. The only new bit is
-tr 0.1 0.1, which sets the target resolution in real-world units. In this case these are degrees of latitude and longitude (as revealed by
gdalinfo). Other common options are meters or feet (which you might run into with some projections in the U.S.)
gdalwarp & the Mercator Projection
gdalwarp -t_srs EPSG:3395 -r lanczos -wo SOURCE_EXTRA=1000 -co COMPRESS=LZW NE1_50M_SR_W_tenth.tif NE1_50M_SR_W_tenth_mercator.tif
Breaking down the code:
gdalwarp invokes the command, while
-t_srs EPSG:3395 sets the target source reference system to EPSG:3395, which is the catalog number for Mercator (I’ll go into other ways to do this in a bit). There are a few ways to find these, spatialreference.org is especially helpful because it displays the description of a map projection in several formats.
-r lanczos defines the resampling method, with a few more options than
gdal_translate. Lanczos is slow but high quality.
-wo SOURCE_EXTRA=1000 is an example of a warp option—advanced parameters that determine how the reprojection is calculated.
SOURCE_EXTRA adds a buffer of pixels around the map as it is reprojected, which helps prevent gaps in the output. Not all reprojections require it, but it doesn’t hurt to add the option to be on the safe side.
-co COMPRESS=LZW works just the same as it does in
NE1_50M_SR_W_tenth.tif NE1_50M_SR_W_tenth_mercator.tif are the input and output filenames, respectively.
Run the command, and you should see a map that looks like this (but slightly larger):
You may notice that the resulting map doesn’t go all the way to 90˚ north or 90˚ south—the Mercator projection becomes infinitely large at the poles, so GDAL is making some guesses about how big the reprojected map should be. It’s part of the magic of GDAL! Which isn’t without its downside. GDAL often behaves unpredictably—in this case, if you reproject the original Natural Earth TIFF instead of the smaller one the resulting Mercator will be a bit taller, stretching from 88.2˚ north to 88.2˚ south. I’ll cover this more in part 3, but add the commands
-te -180 -88.2 180 88.2 -te_srs EPSG:4326 to force a taller map. It’s often good practice to specifiy output extents, even for a global map, since it can save GDAL some guesswork (which can be slow, or in some cases, end badly).
Unfortunately, the Mercator projection has a glaring error at global scale: the further north or south you go the bigger things get. Greenland appears as big as Africa, while in reality it’s only about 7 percent as large. And I won’t even mention Antarctica, which is absolutely huge despite missing its interior.
So instead of making a map that preserves angles, let’s make a map that preserves area. These are especially useful for thematic maps—i.e. maps that display data geographically, not just geography.
The Mollweide Projection
Mollweide is one of many equal-area projections, but has the additional nice property of maintaining straight lines of latitude (also called parallels because they run directly east-west and never meet). We could try to re-reproject the Mercator map to Mollweide, but it’s best to minimize the transformations you apply to a dataset, so let’s work with the rectangular map again. Run:
gdalwarp -t_srs '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -r lanczos -dstalpha -wo SOURCE_EXTRA=1000 -co COMPRESS=LZW NE1_50M_SR_W_tenth.tif NE1_50M_SR_W_tenth_mollweide.tif
(Note: this will generate lots of errors, it’s just GDAL trying to figure out what to do with all the empty space in the corners.)
This is largely the same as the series of commands to make a Mercator map, but I’ve specified the target spatial reference system,
-t_srs , piece by piece instead of relying on an EPSG code. For reasons I don’t understand, some projections are un-loved by the standards bodies, so they don’t have EPSG codes. Mollweide is one of these, so you need to set the parameters manually, using the proj.4 syntax.
In this case
+proj=moll sets the projection to Mollweide, which is followed by a string of settings that define parameters like the center of the map (
+lon_0 (try setting this to 175 if you’re upset that New Zealand always gets short shrift on world maps)).
Rather than describing each property, I’ll just recommend that you find the definitions for projections on the Spatial Reference site. It’s got a nice search function, and defines each projection more than a dozen different ways. On the Mollweide page, for example, you’ll see a list with “proj4” on it—click that link to get the full definition:
+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs. Enclose it in single quotes so it’s parsed correctly.
Finally, I’ve added the option
-dstalpha so the areas that are (literally) off the map are made transparent, instead of given a fill color (black by default).
The Polar Stereographic Projection
Well, no, not perfect—every map is a compromise, remember? What if we wanted to focus on the poles? They’re undefined in Mercator, and shunted off to the edges on most other global maps, including Mollweide. Fortunately, there’s a special map projection designed specifically for the Arctic and Antarctic: polar stereographic.
Some pre-processing is required to create a polar stereographic map using GDAL, which also gives me the opportunity to introduce the concept of a virtual dataset (VRT). A vrt is a text file that functions exactly like a georeferenced data file, but it’s much much smaller. Polar stereographic maps become undefined towards the other pole, so the source data needs to be cropped down to the latitude that will become the edge of the map. I’ll again use
gdal_translate, but with some new options:
gdal_translate -of VRT -projwin -180 -60 180 -90 NE1_50M_SR_W_tenth.tif NE1_50M_SR_W_SH60.vrt
Instead of letting
gdal_translate write the default GeoTIFF, I’ve used
-of to specify the output format as a VRT. I’ve then cropped out Antarctica with
-projwin specifying the boundary in the following format: upper left x, upper left y, lower right x, lower right y (but without the commas). In the Southern Hemisphere upper left longitude is -180˚, upper left latitude is -60˚, lower right longitude is 180˚, and lower right latitude is -90˚ (the South Pole).
With a nice compact (and quick to make) VRT, here’s the command to make a map of Antarctica:
gdalwarp -t_srs EPSG:3976 -ts 7200 0 -r near -dstalpha -wo SOURCE_EXTRA=1000 -co COMPRESS=LZW NE1_50M_SR_W_SH60.vrt NE1_50M_SR_W_sh60_polarstereo.tif
It starts like the command for Mercator, specifying the target spatial reference system with an EPSG code:
EPSG:3976. Then it gets (very slightly) weird. I’ve forced the target size to be 7,200 pixels with
-ts 7200 0 and then set resampling to nearest neighbor with
-r near. This gets around a nasty issue that was blurring pixels along the +/−180˚ line, at the expense of making the map look pixellated (at least if you zoom in). The rest mirrors the previous commands.
I set the target size to be so much larger than the resolution needed for this post (only 1,400 pixels across) so that I could reduce the resolution with a resampling method that averages pixels—bilinear. Here’s a
gdal_translate command to shrink the image and save it as a PNG, a nice lossless compression method (like GIF) that also retains full color (like JPEG) and displays on the Web.
gdal_translate -of PNG -outsize 1400 0 -r bilinear NE1_50M_SR_W_sh60_polarstereo.tif NE1_50M_SR_W_sh60_polarstereo_1400.png
A nice map of Antarctica extending from the South Pole to 60˚ north. I’ll leave it as an exercise for the reader to make a similar map of the Northern Hemisphere.
One final note: although GDAL is very powerful and very flexible, it can’t do everything in mapping. Some projections, even some of my favorites like Hammer-Aitoff, Winkel-tripel, and Natural Earth, are not fully supported. Before trying out your favorite projection check the Spatial Reference database and if the proj.4 entry is blank, it probably won’t work. The good news is support for some new map projections is underway.
- NASA’s DSCOVR: EPIC page with real-time views of the sunlit Earth.
- JSC’s Earth Observation Lab for a closer view.
- The David Rumsey Map Collection for a historical perspective.
- Natural Earth for comprehensive raster & vector data.
- Axis Maps’ Thematic Cartography Guide.
- Spatial Reference.
- Hands on with GDAL/OGR: a gentle introduction to command line GIS.
- Even Rouault’s GDAL workshop.
- XKCD on Map Projections.
- A Gentle Introduction to GDAL
- Map Projections & gdalwarp (you are here)
- Geodesy & Local Map Projections