4. Tutorial 1: simple raster conversion with GDAL

For this tutorial, you will need to have the GDAL raster utilities. You can download these from the GDAL website.

4.1. First, get a DEM (digital elevation model)

You can get DEMs from all over the place, here are some common data distribution centres:

Much more data can be found via search engines.

4.1.1. The DEM used in some of the examples here

The DEM used in this example can be downloaded here.

4.2. Now start a terminal window

Instructions for starting a terminal window are here: Getting onto our servers using NX.

You will need to navigate to the place you’ve got the DEM using commands like:

  • cd and cd .. to change directories.
  • ls to see what files are in the current directory.
  • pwd to see what folder you are currently in.

If you’ve correctly navigated to the folder and run ls and pwd commands you will get something a bit like this:

smudd@burn Scotland $ pwd
smudd@burn Scotland $ ls
smudd@burn Scotland $

4.3. Seeing what the DEM contains using gdalinfo

The GDAL tool gdalinfo is used to see if GDAL can actually recognise your data a a geographic dataset. You can try it out on the test DEM:

smudd@burn Scotland $ gdalinfo WhiteadderDEM.tif

will give you all sorts of information:

Driver: GTiff/GeoTIFF
Files: WhiteadderDEM.tif
Size is 486, 645
Coordinate System is:
PROJCS["Transverse Mercator",
    GEOGCS["OSGB 1936",

and so on. Basically if this stuff gets spit out, then GDAL thinks it is geographic raster and you probably can read it using a GIS like ArcMap or QGIS.

4.4. Converting between formats using GDAL

In this tutorial we will convert the format of the DEM. The reason why is because the University of Edinburgh topographic analysis software takes specific data formats. All three can be read by GIS software:

  • asc An ascii-based data format. Basically you can open up these files in a text editor and look at the data but the files are HUGE so we almost never use this.
  • flt This is a binary data format so you can’t just look at the data in a text editor but the files are much smaller. This file format has the disadvantage that is doesn’t contain georeferencing information. Georeferencing is just a fancy name for all the data that tells a GIS where in the world your data is.
  • bil The third option is something called an ENVI bil raster which has the extension bil. This is the preferred option since it is both small (versus the ascii format) and retains georeferencing.
  • Both the flt and bil formats include a separate hdr file that contains some information about the DEM.

We can convert files using the GDAL tool gdal_translate.

gdal_translate recognises many file formats (see list here), but for LSDTopoTools you want either:

  • The ESRI .hdr labelled format, which is denoted with EHdr.
  • The ENVI .hdr labelled format, which is denoted with ENVI. ENVI files are preferred since they work better with GDAL and retain georeferencing.

To set the file format you use the -of flag, an example would be:

$ gdal_translate -of ENVI WhiteadderDEM.tif WhiteadderDEM.bil

Okay, lets see the result:

smudd@burn Scotland $ gdal_translate -of ENVI WhiteadderDEM.tif WhiteadderDEM.bil
Input file size is 486, 645
0...10...20...30...40...50...60...70...80...90...100 - done.
smudd@burn Scotland $ gdalinfo WhiteadderDEM.bil
Driver: ENVI/ENVI .hdr Labelled
Files: WhiteadderDEM.bil
Size is 486, 645
Coordinate System is:
    GEOGCS["GCS_OSGB 1936",

and a bunch more georeferencing information.


4.5. Clipping rasters using GDAL

You might also want to clip your raster to a smaller area. This can sometimes take ages on GUI-based GISs. An alternative is to use gdalwarp for clipping:

$ gdalwarp -te <x_min> <y_min> <x_max> <y_max> input.bil clipped_output.bil

Since this is a gdalwarp operation, you can add all the bells and whistles to this, such as:

  • changing the coordinate system,
  • resampling the DEM,
  • changing the file format.

More documentation on gdalwarp can be found on the official site <http://www.gdal.org/gdalwarp.html>_ or here at the LSDTopoTools documentation page.

The main thing to note about the -te operation is that the clip will be in the coordinates of the source raster (input.bil). You can look at the extent of the raster using gdalinfo.

4.6. Merging rasters using GDAL

This is called mosaic to new raster in ArcMap. In ArcMap it is VERY slow. If you have large DEMs you can waste huge amounts of time trying to do this in ArcMap. If you are going to merge rasters, trust me, GDAL is MUCH better.

You merge rasters in GDAL using gdal_merge.py.

Usage looks something like this:

gdal_merge.py -o out.tif in1.tif in2.tif

The -o indicates the outfile. Don’t forget!

If you want to merge to an ENVI file (the recomended file format) you could do this:

gdal_merge.py -o out.bil -of ENVI in1.tif in2.tif

You use the -of flag to denote the output format (in this case an ENVI file).

If you want to merge and at the same time change the pixel size, you do this:

gdal_merge.py -o out.bil -of ENVI -ps 30 30 in1.tif in2.tif

Here the -ps flag is for pixel size. Assuming this is a projected coordinate system (see below) in metres, you will get 30x30m pixels. Be careful with this: the GDAL documentation doesn’t say what sort of interpolation method is used when resampling, so if you want to change the pixel size, use gdalwarp.

4.7. Changing raster projections with gdalwarp

The preferred coordinate system is WGS84 UTM coordinates. For convert to this coordinate system you use gdalwarp. The coordinate system of the source raster can be detected by gdal, so you use the flag -t_srs to assign the target coordinate system. Details about the target coordinate system are in quotes, you want:

'+proj=utm +zone=XX +datum=WGS84'

where XX is the UTM zone. Here is a map of UTM zones. For example, if you want zone 44 (where the headwaters of the Ganges are), you would use:

'+proj=utm +zone=XX +datum=WGS84'

Put this together with a source and target filename:

$ gdalwarp -t_srs '+proj=utm +zone=XX +datum=WGS84' source.ext target.ext

so one example would be:

$ gdalwarp -t_srs '+proj=utm +zone=44 +datum=WGS84' diff0715_0612_clip.tif diff0715_0612_clip_UTM44.tif

There are several other flags that could be quite handy (for a complete list see here).

  • -of format: This sets the format to the selected format. This means you can skip the step of changing formats with gdal_translate. We will repeat this later but the formats for LSDTopoTools are

    • AAIGrid for ASCII files. These files are huge so try not to use them.
    • EHdr for ESRI float files. This used to be the only binary option but GDAL seems to struggle with it and it doesn’t retain georeferencing.
    • ENVI for ENVI rasters. This is the preferred format. GDAL deals with these files well and they retain georeferencing. We use the extension bil with these files.

    So, for example, you could output the file as:

    $ gdalwarp -t_srs '+proj=utm +zone=44 +datum=WGS84' -of ENVI diff0715_0612_clip.tif diff0715_0612_clip_UTM44.bil
  • -tr xres yres: This sets the x and y resolution of the output DEM. It uses nearest neighbour resampling by default. So say you wanted to resample to 4 metres:

    $ gdalwarp -t_srs '+proj=utm +zone=44 +datum=WGS84' -tr 4 4 diff0715_0612_clip.tif diff0715_0612_clip_UTM44_rs4.tif

    IMPORTANT: LSDRasters assume square cells so you need both x any y distances to be the same

  • -r resampling_method: This allows you to select the resampling method. The options are:

    • near: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
    • bilinear: bilinear resampling.
    • cubic: cubic resampling.
    • cubicspline: cubic spline resampling.
    • lanczos: Lanczos windowed sinc resampling.
    • average: average resampling, computes the average of all non-NODATA contributing pixels. (GDAL >= 1.10.0)
    • mode: mode resampling, selects the value which appears most often of all the sampled points. (GDAL >= 1.10.0)

    So for example you could do a cubic resampling with:

    $ gdalwarp -t_srs '+proj=utm +zone=44 +datum=WGS84' -tr 4 4 -r cubic diff0715_0612_clip.tif diff0715_0612_clip_UTM44_rs4.tif
  • -te <x_min> <y_min> <x_max> <y_max>: this clips the raster. You can see more about this above in Clip-using-gdal.