Astronaut photograph ISS050-E-66060 courtesy NASA Earth Observation Lab.

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, gdalinfo, and gdal_translate.)

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.

The first photograph showing a curve on Earth’s horizon, taken on November 11, 1935. Courtesy Jerry Bryant, Lawrence County Historical Society.

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.

Large parts of California and Nevada are visible in this photograph taken from the Internation Space Station on October 19, 2014. Astronaut photograph ISS041-E-081461 courtesy NASA Earth Observation Lab.

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.)

The Western Hemisphere on April 7, 2017. From DSCOVR’s EPIC instrument.

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.

The earliest known globe was crafted in 1492 — the same year Columbus sailed to America (but before he came back—note the lack of the Americas). To construct a globe this printed map would have been cut out and folded into a sphere, with the poles pasted on. Composite: Globe Gores 1–4. Martin Behaim’s Erdapfel, 1492. Courtesy David Rumsey Map Collection.

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.

Mercator, Mollweide, and Southern Hemisphere polar stereographic maps. Made with Natural Earth.

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.

Natural Earth I with Shaded Relief and Water. This is in a very common projection with way too many names: Plate Carrée, equirectangular, equidistant cylindrical, simple cylindrical, rectangular, lat-lon, geographic projection, WGS 84, or EPSG:4326. It’s simply an even grid of latitude and longitude, centered at 0˚north 0˚south. Made with Natural Earth.

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

Now let’s convert this map from Plate Carrée/rectangular/lat-lon/etc. to Mercator, using another GDAL utility, gdalwarp:

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 gdal_translate, and 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):

An (almost) global map using the Mercator projection. Since the Mercator projection stretches to infinity at the poles, the highest northern and southern latitudes are automatically clipped via mysterious GDAL guesstimates. Made with Natural Earth.

Success!

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 Mollweide projection, which is equal-area and has the bonus feature of maintaining straight lines of latitude. Made with Natural Earth.

Perfect!

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 polar stereographic projection, centered on the South Pole and extending to 60˚ south. Notice how razor-sharp the edge is—that’s the result of the render large then downsize technique, combined with lossless compression. Made with Natural Earth.

Huzzah!

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.

Phew. Next up: geodesy, local map projections & gdal_merge.py.

Further Reading

• 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.
• GDAL.org.
• Hands on with GDAL/OGR: a gentle introduction to command line GIS.
• Even Rouault’s GDAL workshop.
• XKCD on Map Projections.

1. A Gentle Introduction to GDAL
2. Map Projections & gdalwarp (you are here)
3. Geodesy & Local Map Projections
4. Working with Satellite Data
5. Shaded Relief
6. Visualizing Data
7. Transforming Data
8. Reading Scientific Data Formats

Thanks to Frank Warmerdam, Kelsey Jordahl, and Joe Kington for technical support and editing.