More on Processing UK Environment Agency LIDAR Data

Since I found that my last blog post on displaying Environment Agency LIDAR DEMs has become the most ever viewed blog post that I have written (popular subject apparently), I've been thinking of writing a followup, having learned a few new things. One of the main problems with dealing with all this LIDAR data is speed. First, getting the data to a state that is useful takes a lot of processing, which can be solved by automating the process using some python scripts I wrote. Second, the draw speed on the screen in QGIS can be solved using Virtual Raster Tables and Pyramids. This tutorial will assume that you are using the OSGeo4W package on windows, with QGIS, python  and OSGeo4W Shell options installed, but much of it may be transferable to other setups.

I'll start by throwing some python code at you for processing the data, and then explain a bit about what it does and how to run it. Put this code in a file somewhere called demimport.py


#!/usr/bin/env python

import sys
import zipfile
import os

impdir='e:\\data\\gis\\dem\\import'
expdir='e:\\data\\gis\\dem'

def unzipfile(filename,exportto):
    with zipfile.ZipFile(filename, "r") as z:
        z.extractall(exportto)

def main( argv=None ):
    # Part 1, unzip the data and merge the contents into one file, deleting the intermediate files on completion
    for file in os.listdir(impdir):
        if file.endswith(".zip"):
            # Work out the filenames we are using
            print "File: ", file
            stubname = file[:-4]
            print "Stubname: ", stubname
            outdir = expdir + "\\" + stubname
            print "outdir: ", outdir
            impfile = impdir + "\\" + file
            print "impfile: ", impfile
            demname = outdir + "\\" + stubname + ".asc"
            print "demname: ", demname

            # First, deal with the zip file
            if not os.path.exists(outdir):
                os.mkdir(outdir)
            unzipfile(impfile,outdir)
            os.remove(impfile)

            # Now get a list of files in the created directory and construct the script to merge them
            script = "gdal_merge.bat -n -9999 -a_nodata -9999 -of GTiff -o " + demname
            for file2 in os.listdir(outdir):
                if file2.endswith(".asc") and not file2.startswith("LIDAR-DTM"):
                    script += " " + outdir + "\\" + file2
            print "script: ", script
            os.system(script)
            for file2 in os.listdir(outdir):
                if file2.endswith(".asc") and not file2.startswith("LIDAR-DTM"):
                    os.remove(outdir + "\\" + file2)

    # Part 2, Create hillshade files from the files created in part 1
    for f in os.listdir(expdir):
        if os.path.isdir(os.path.join(expdir, f)) and f.startswith("LIDAR-DTM"):
            # Work out the filenames we are using
            stubname = expdir + "\\" + f + "\\" + f
            print "Stubname: ", stubname
            demname = stubname + ".asc"
            print "demname: ", demname
            hsname = stubname + "-HS.asc"
            print "hsname: ", hsname
            hs2name = stubname + "-HS2.asc"
            print "hs2name: ", hs2name

            if os.path.exists(demname):
                if not os.path.exists(hsname):
                    # Now create the hillshade
                    script = "gdaldem hillshade " + demname + " " + hsname + " -z 1.0 -s 1.0 -az 315.0 -alt 45.0 -compute_edges -of GTiff"
                    os.system(script)
                if not os.path.exists(hs2name):
                    # Now create the second hillshade
                    script = "gdaldem hillshade " + demname + " " + hs2name + " -z 1.0 -s 1.0 -az 45.0 -alt 45.0 -compute_edges -of GTiff"
                    os.system(script)

if __name__ == '__main__':
    sys.exit(main())

You will see near the top of this script two directories called impdir and expdir. impdir is the directory where you dump the zipfiles you download from the environment agency website. expdir is the directory where the script will output the resulting data. Change these to whatever directory structure you wish to use on your computer. The script will  merge the contents of each zip file into a single file and then create two different hillshades from that data. More on the hillshades later. It automates most of what I described in my last blog post, all apart from setting the style data in QGIS. You can run the script by opening OSGeo4W Shell, changing directory to where you have saved the python file and typing @python demimport.py.

The next bit is about speeding up the data display in QGIS itself, using the aforementioned Virtual Raster Tables (hereafter, VRTs) and Pyramids. First, VRTs. Imagine you have 100 LIDAR tiles active for display in QGIS. Each time you zoom or move the display, it has to read them all to see which it can display in the area you are viewing at the time, which obviously means a lot of slow reading from your hard disk (assuming you still use those). A VRT acts as an index file for your 100 LIDAR tiles, so QGIS only has to look in one place to find out what it needs to display, and then only has to open the tiles that it requires to fill the area you are viewing, so everything displays much quicker. Pyramids speed things up in a different way. Imagine you are quite zoomed out, looking at a wide area of LIDAR data. QGIS would normally have to read the whole of each file and reduce it in size to fit into the small area that the tile would be displayed in. Creating a pyramid does this process ahead of time, so it takes the original image, compresses ito to a quarter of its size, then does that again and again and saves all that in a pyramid file, so when you are zoomed out, rather than reading the entirity of the original data, it will read the pre-compressed data suitable to your zoom level, reducing the amount it has to read from your hard disk and speeding up the display process. Those are the concepts behind it, now for the code that actually does it. Put this code in a file called makevrt.py

#!/usr/bin/env python

import sys
import os



def main( argv=None ):
    # definitions for the vrt we want to create, this changes
    resolution = "2M"
    expdir = 'e:\\data\\gis\\dem'
    inputareas = ["ST","SY","SS","SX"]

    # File and directory names, this doesn't change
    ldtm = "LIDAR-DTM-"
    ext = ".vrt"
    vrtnamestub = "VRT" + "\\" + ldtm + resolution
    types = [ "", "-HS", "-HS2" ]    listname = expdir + "\\filelist.txt"

    for type in types :
        for ia in inputareas:
            vrtname = expdir + "\\" + vrtnamestub + "-" + ia + type + ext
            pyramidname = vrtname + ".ovr"

            if not os.path.exists(vrtname):
                script = "gdalbuildvrt -input_file_list " + listname + " " + vrtname

                base = ldtm + resolution + "-" + ia
                base2 = expdir + "\\" + base
                fp = open(listname,"w")
                for dir in os.listdir(expdir):
                    if dir.startswith(base):
                        fp.write(expdir + "\\" + dir + "\\" + dir + type + ".asc\n")
                fp.close()

                print "script: " + script
                os.system(script)

            if os.path.exists(vrtname) and not os.path.exists(pyramidname):       
                print "processing: " + vrtname
                # Now create the pyramids
                script = "gdaladdo -r CUBIC -ro " + vrtname + " 2 4 8 16 32 64"
                os.system(script)


if __name__ == '__main__':
    sys.exit(main())

You can run the script by opening OSGeo4W Shell, changing directory to where you have saved the python file and typing @python makevrt.py. Again, there are some changes you will need to make at the top of the file. resolution is the type of data you are dealing with. You can download 2M, 1M, 50CM or 25CM data from the Environment Agency website, this script will only do one at a time. expdir should be the same as expdir from the first script. You will also need to create a subdirectory off of that called 'VRT', which is where this script creates its new files. inputareas is a list of Ordnance Survey Grid Letters that the script will generate VRTs for. It will generate one set of files for each grid letter, you have to tell it which ones to do. After that, you can load the newly created VRT files into QGIS using the Add Raster Layer button and style them as per my original blog post.

Now back to those two different hillshades I mentioned. Why two different hillshades? If you imagine a hillshade a shining a light from a particular angle to create highlights and shadows, a linear feature running in the same direction as the direction of the light will not show up very well, so I've created a second hillshade with the light coming in at a 90 degree angle to the first hillshade. You can see the difference below.

 

 First hillshade. Click for larger image.

 Second hillshade. Click for larger image.

If you click on the images and look at the highlighted feature, a possible new (unconfirmed) Roman road on the Isle of Wight, you will see that it is much clearer in the second image compared to the first. My script generated the hillshades with light coming in from the north-east in the first image, which is parallel to the road feature, and from the north-west in the second image, which is perpendicular to the road feature, showing it up better.