================================================================ 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. First, get a DEM (digital elevation model) ============================================ You can get DEMs from all over the place, here are some common data distribution centres: * For LiDAR, two good sources are `opentopography `_ and the `U.S. Interagency Elevation Inventory `_. Other countries are not so progressive about releasing data but you could rummage around the `links here `_. * For 10m data of the United States you can go to the `national map viewer `_. To get ordinance survey data at 10m resolution you should use `Digimap `_ note: requires institutional subscription). * ASTER 30m data can be obtained from `NASA's reverb site `_. * One site that has filled and corrected `90-m DEMs from the SRTM mission `_. Much more data can be found via search engines. The DEM used in some of the examples here ------------------------------------------- The DEM used in this example can be downloaded `here `_. Now start a terminal window ============================================ Instructions for starting a terminal window are here: :ref:`getting-onto-servers`. 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 /home/my_sensibly_named_folder/Scotland_Data smudd@burn Scotland $ ls WhiteadderDEM.tif smudd@burn Scotland $ 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. .. _change-formats-GDAL: 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 WhiteadderDEM.bil.aux.xml WhiteadderDEM.hdr Size is 486, 645 Coordinate System is: PROJCS["Transverse_Mercator", GEOGCS["GCS_OSGB 1936", and a bunch more georeferencing information. .._Clip-using-gdal: 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 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 _` 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``. 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``. 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`` ` `: this clips the raster. You can see more about this above in :ref:`Clip-using-gdal`.