Sunset at the Dollhouse. Maze District, Canyonlands National Park. Photo courtesy National Park Service.

A Gentle Introduction to GDAL, Part 1

Published in

Planet Stories

9 min readApr 4, 2017

In 1885 the United States Geological Survey published their first map of the Canyonlands, a remote, rugged, high, dry, and otherwise inaccessible chunk of desert in southeastern Utah.

USGS 1:250000-scale Quadrangle for La Sal, UT 1885. Via the National Geologic Map Database.

These maps were based on the “triangulation and topography” recorded during a pair of expeditions lead by Major John Wesley Powell, who later become a director of the Survey.

In a roadless landscape heavily sculpted by erosion, the easiest way through was via the rivers that wound their way through steep-walled canyons, so Powell and his men rafted down the Green and Colorado Rivers (then called the Grand River). Over the course of two journeys in 1869 and 1871, the Powell expeditions carefully mapped the landscape as they went.

Unfortunately, they were largely stuck at the bottom of those canyons: a hike to the top to get a better view was usually arduous and often impossible. Once at the rim of the canyons the men were met with inhospitable slickrock, carved into fantastic shapes, which further inhibited travel and blocked sight-lines.

On the top of the country immediately at the junction of the Grand and Green Rivers, westside. Utah. Photograph taken by the second Powell expedition, September 17, 1871. Photo courtesy USGS.

The resulting maps were largely accurate with respect to the courses of major rivers and the locations of high peaks, but were lacking in fine detail. Largely absent were the slot canyons, grabens, arroyos, arches, pinnacles, salt domes and other geological curiosities that make the area so fascinating today.

Not that any blame lies on Major Powell or his crew: the area is so inaccessible and difficult to traverse that the USGS didn’t substantially update its maps for over 65 years—new printings were issued in 1896, 1901, 1909, 1918, and 1923 with only minor typographic changes. Large-scale maps of the Canyonlands weren’t published until the 1950s—among the last parts of the continental United States to be charted in detail.

Comparison of the USGS 1:250,000-scale topographic maps of the confluence of the Green and Colorado Rivers drawn in 1885 (1923 printing, left) and 1959 (right). From the USGS National Geological Map Database.

Two factors led to these radical improvements after decades of stasis—advances in photogrammetry which enabled map-makers to precisely determine topography from overlapping aerial photographs, and a post-war uranium boom that inspired prospectors to crawl over ever inch of the convoluted landscape.

Aerial photograph of the Canyonlands taken by the USGS in 1953. The Green River merges with the Colorado River in the center of the frame. Photograph from the USGS Earth Explorer.

It’s been another 65 years since the terrain of the Canyonlands was finally mapped, and imagery from above is now ubiquitous. Once-remote parts of the globe are now pictured by aircraft, satellites and drones; using visible light (red, green, blue and the other colors of the rainbow), infrared light, thermal radiation, LIDAR, and radar from a scale of centimeters to kilometers. Weather satellites capture images of an entire hemisphere every 15 minutes, and Planet (my company) will soon be imaging the entire surface of the Earth (well, the cloud-free bits at least) at about 4 meters per pixel (the size of a small building) every day.

Accompanying this explosion of raster (pixel-based) data was growth in the amount and availability of vector data—everything from county boundaries to geological maps to GIS locations attached to photographs on flickr. This bewildering array of data is accompanied by a bewildering array of data types: file formats, scales, map projections, datums, etc., etc., etc.

Locals and Tourists, Boston. Blue pictures are by locals. Red pictures are by tourists. Yellow pictures might be by either. ©2010 Eric Fisher, CC BY-SA 2.0.

Why GDAL?

Given that these data are (usually) much more useful combined than separate, how does one integrate disparate data into interesting and novel composites?

Broadly speaking, the answer is a GIS—geographic information system (perhaps the first and last time I link directly to a Wikipedia article). “A system designed to capture, store, manipulate, analyze, manage, and present spatial or geographic data.”

The bad news is that the most pervasive GIS software is expensive, difficult to learn, and won’t even run on my operating system of choice. The good news is that there’s an open-source alternative—GDAL—that’s free, broadly supported, constantly updated, and runs on almost anything. Excepts it’s also difficult to learn, especially if you’re terrified of the command line, like me.

GDAL stands for “Geospatial Data Abstraction Library”:

A translator library for raster and vector geospatial data formats that is released under an X/MIT style Open Source license by the Open Source Geospatial Foundation. As a library, it presents a single raster abstract data model and single vector abstract data model to the calling application for all supported formats. It also comes with a variety of useful command line utilities for data translation and processing.

Even the description is complicated. Put another way, GDAL (pronounced ‘gee-dal’ not ‘goodle’ (too close to Google!)) is a collection of software that helps with the translation of data from different file formats, data types, and map projections. It’s used in free software like QGIS and GRASS, and even some commercial applications like ESRI ArcGIS and Avenza Geographic Imager/MAPublisher.

Installing GDAL

As I already mentioned GDAL runs on the command line, but is only a little bit harder to install than the typical commercial app.

OSX

On a Mac, I’ve found the easiest approach is to download the GDAL 2.1 Complete disk image from KyngChaos (I know it sounds like a mid-80s hacker collective distributing bootleg C64 games, but the site’s legit). This installer relies on the version of Python included with OS X, in my case 2.7.12, not 3.x.

On the disk image are two installers: one for NumPy (run this first) and one for GDAL Complete (run this second). You’ll have to open the files with a control-click since the installers are from an “unknown developer.” (Don’t worry, it’s still ok. (I think.)) After one last step you’ll be good to go: cut and paste the following two lines (borrowed from MapBox) into the Terminal (they’ll allow you to run GDAL commands just by typing them, without specifying a path. bash_profile is a hidden file, so don’t be surprised if you can’t find it.)

echo 'export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH' >> ~/.bash_profilesource ~/.bash_profile

Windows

Follow these instructions to install both Python 2.7 and GDAL 2.1 from the UCLA Institute for Digital Reseach and Education. (Be careful to set the paths correctly for a 32- or 64-bit Windows install.)

Linux

If you’re running Linux you know more about this than I do—you’re on your own.

An alternative which works on all platforms is to use the the Conda package manager. If you’re already running Conda you can install GDAL with:

conda install -c conda-forge gdal=2.1.3

First Steps Using GDAL

To make sure GDAL is installed, open up a command line window and type:

gdalinfo --version

If you see something like:

GDAL 2.1.2, released 2016/10/24

You’re good! (Hmm, it looks like I’m a tiny bit out of date.)

Now, let’s download a map to work with—a shaded relief of Canyonlands, cropped from the U.S. National Park Service map (which has a fantastic cartography center, BTW).

Shaded relief map of the Canyonlands. Courtesy National Park Service Harpers Ferry Center.

gdalinfo

Now that you’ve got an image, navigate to the folder you downloaded to and enter this into the command line:

gdalinfo CANYrelief1-geo.tif -mm

“gdalinfo” runs the eponymous utility program (one of many included with GDAL), “CANYrelief1-geo.tif” is the name of the file you just downloaded, and “-mm” is a command that calculates and displays additional information. After hitting return, you should see a block of text, beginning with:

Driver: GTiff/GeoTIFF

This indicates that the file is a GeoTIFF, which is a special type of TIFF that stores the information necessary to place each pixel in the image on the surface of the Earth. It’s an incredibly flexible and useful format, and is increasingly by adopted to store more types of data than just images.

Reading further, the next lines of text are:

Files: CANYrelief1-geo.tifSize is 2800, 2800

These simply indicate the file name, and the size (in pixels) of the image. Here’s the really important bit:

Coordinate System is:PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]Origin = (-12249462.599999999627471,4629559.794860946945846)Pixel Size = (13.284000000000001,-13.285397060378999)

This is the information that places this image in the Canyonlands, and specifies the location of each pixel. I’ll get into map projections in my next post, but if you want to learn more the USGS has a nice primer (PDF)—or just read XKCD. The last entry is the pixel size in meters—defined by the “unit” entry further up.

Next up is some generic metadata for the file:

Metadata:AREA_OR_POINT=AreaTIFFTAG_DATETIME=2017:04:01 20:24:57TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)TIFFTAG_SOFTWARE=Adobe Photoshop CC (Macintosh)TIFFTAG_XRESOLUTION=72TIFFTAG_YRESOLUTION=72Image Structure Metadata:COMPRESSION=LZWINTERLEAVE=PIXEL

This tells me that I made the image on April Fool’s Day, I used Photoshop to crop and save the file (with a plugin that keeps the geographic info intact), the resolution (if printed) is 72 dots per inch, and the image is saved with LZW, a type of lossless compression.

Corner Coordinates:Upper Left  (-12249462.600, 4629559.795) (110d 2'19.66"W, 38d21'14.51"N)
Lower Left  (-12249462.600, 4592360.683) (110d 2'19.66"W, 38d 5'29.42"N)
Upper Right (-12212267.400, 4629559.795) (109d42'16.79"W, 38d21'14.51"N)
Lower Right (-12212267.400, 4592360.683) (109d42'16.79"W, 38d 5'29.42"N)
Center      (-12230865.000, 4610960.239) (109d52'18.23"W, 38d13'22.39"N)Band 1 Block=2800x31 Type=Byte, ColorInterp=RedComputed Min/Max=149.000,255.000Band 2 Block=2800x31 Type=Byte, ColorInterp=GreenComputed Min/Max=133.000,255.000Band 3 Block=2800x31 Type=Byte, ColorInterp=BlueComputed Min/Max=100.000,255.000

The last block of text displayed by gdalinfo shows the corner points of the image in two different units (meters and minutes, degrees, seconds) and finally some information about each band (also called a channel )in the image. Each channel is byte (8-bit) format, there are three bands (red, green, and blue, respectively), and the minimum and maximum values in each band.

That last detail was calculated and displayed because I added “-mm” (min/max) to the command. If you throw in “-stats” you’ll get additional statistics like mean, median, and standard deviation. But be careful, some software (QGIS) may write estimates of these values into the file header and gdalinfo will return the wrong values if you just use “stats”, whereas “-mm” will force statistics to be calculated on the full dataset. Slow but thorough.

gdal_translate

Now that we have some information about the file, let’s do something useful—resizing the image to a web-friendly 1,400-pixels-wide and saving it as a JPEG—with “gdal_translate”. Enter these commands into the terminal:

gdal_translate -of JPEG -co QUALITY=70 -co PROGRESSIVE=ON -outsize 1400 0 -r bilinear CANYrelief1-geo.tif CANYrelief1.jpg

You’ll see two lines of text as a result, and have a brand-new file, CANYrelief1.jpg, in the same directory as the TIFF. I specified the file type, size, etc. with the following commands:

-of JPEG -co QUALITY=70 -co PROGRESSIVE=ON

“-of” sets the output format, in this case JPEG. “-co QUALITY=90” and “-co PROGRESSIVE=ON” are creation options, a series of commands specific to each file type, in this case writing a JPEG with a quality of 90 out of 100 that will load progressively.

-outsize 1400 0 -r bilinear

These two commands indicate set the output size and interpolation method of the resizing. “-outsize” is set in pixels—x (horizontal) first and y (vertical) second. In this case I set y to zero, so the image kept it’s original aspect ratio. I would have gotten the same results with “-outsize 1400 1400”. “-r” specifies the resampling method, of which GDAL offers a generous selection (nearest, bilinear, cubic, cubicspline, lanczos, average, mode). “Nearest” is the default, but it leaves visible stairsteps so you’ll almost always want to change it. I usually use “bilinear” for satellite imagery since it’s quick and I sharpen as an additional step.

Finally:

CANYrelief1-geo.tif CANYrelief1.jpg

These are just the input file name and the output file name, respectively (but be careful, gdal_translate will happily overwrite files with no warning).

gdalinfo and gdal_translate are two of the more straightforward utilities included with GDAL. And, to be honest, there are a thousand ways to convert a TIFF to a JPEG—but not nearly so many are also able to read the headers in a geotiff, translate obscure data formats (if you do data visualization long enough, you’ll run into some special files), or transform a map from one projection to another. Next up is Part 2: Map Projections & gdalwarp, which begins to unlock the power of geospatial data, but requires some familiarity with map projections (spoiler: the Earth isn’t flat).