29 October 2016

Latest Results: Chichester

Following on from last year's very successful radar survey in Chichester, I went back for another week of the same, this time around the area of the Cathedral in the south-west quadrant of the city. Chichester & District Archaeological Society had already found a lot in the area using earth resistance and excavation, so the radar didn't show a lot that was new, just in slightly better definition. There wasn't a lot around the cathedral itself, so the area has probably had any Roman remains there thoroughly removed, but there was around the Deanery and Bishop's Palace. The garden of the Deanery contains the old medieval (or post-medieval) Deanery, and the the area in front of the Bishop's Palace contains a Roman building and the medieval hall that was the old Bishop's palace.

The old Deanery in the garden of the current Deanery
Click for larger image

Roman building (green) and medieval Bishop's Palace (orange)
Click for larger image

I also went back to Priory Park to look again at the third Roman building (in green, within the lighter survey area) found near the cricket pitch. It has suffered greatly from robbing, not least by the Saxons, who seem to have used some of the stone in their sunken floor buildings, two of which appear (in purple) in the higher resolution re-survey of this area.

Priory Park survey. Click for larger image.

I'll be giving a talk on the results from both years for CDAS on the 22nd of February 2017 at 7:30pm in the cinema of the New Park Centre, New Park Road, Chichester.


04 October 2016

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.







14 August 2016

Roman Roads conference in Portsmouth

Many of you who know me know that my personal research project, which I have been working on for the last few years is Roman Roads and roadside settlements. I've been asked to speak on the subject at a conference on the subject of Roman roads in memory of Ivan Margary, the most famous of all Roman road researchers in Britain. I'll be doing two talks (that's a full Bonsall for you geophysicists out there), the first on the subject of Roman roads research in the county of Sussex, Ivan's original stomping ground, and the second talk on the use of geophysics for finding Roman roads. The conference is a two day conference on the 3rd and 4th of September 2016 at Portsmouth University. The website for the organisers of the conference is here.


07 April 2016

Snuffler version 1.2: The radar experiment

It's time for a new version of Snuffler. I've upgraded from visual studio 2005 to 2015 and changed drawing of custom graphics in dialog boxes, e.g. in the Display Settings dialog, so let me know if you get any strange new problems. If you get an error starting up the new version of Snuffler, you will probably need to run the redist exe that comes in the same zip file.

The main feature added this time is functionality for the display of radar data. Not much mind you, very little in fact. It all started back when I got my radar. I also got a copy of ReflexW, and I'm happy enough using that, but having written Snuffler, I found I was not completely happy with the way radar data is traditionally displayed, and wanted to have a go at writing a display process to see if I could do better. Apart from basic training on how to use ReflexW, I have no formal training in GPR processing theory. Thousands of papers on the subject have passed me by, and I wouldn't really know where to start looking. Some of what I'm doing here is re-inventing the wheel, but I still wanted to have a go at it, despite going in blind. The new functionality I have added to Snuffler is not intended to be a useful GPR processing system, it is very far from that. It is intended to scratch an itch I had regarding GPR display methods and learn something for myself along the way.

For those of you who are not completely familiar with how radar data is displayed, here is a bit of history to bring you up to speed. Apologies if you already know this stuff. Click on any of the images to show a larger version. Unlike for earth resistance and magnetometry, which give a single reading at any point, radar will give you a vertical trace going through the ground. The resulting wave recorded at a single point in horizontal space is known as an A-Scan, and an example is shown below.

A-Scan

The top of this wave is at the ground surface and goes down into the ground. The dotted red line is an amplitude (signal strength) of zero, so as the wave extends further from that line, the more signal is being reflected and the greater the strength of signal. You tend to get a greater reflection at the surface, as less and less energy penetrates further into the ground as you go deeper. So that's a single point on the ground. If you record a line of these and string them together, this is known as a B-Scan. An example of which is shown below.


 B-Scan, historical display method

This is how B-Scans used to be displayed in the early days of computing, when everything was monochrome. The A-Scans where arranged in a line, and the area under the curve of the wave to the right of the zero line was filled with black to create contrast for high amplitude areas. This display method is the equivalent of the dot-density plot for earth resistance, when technology was similarly limited to monochrome, and thankfully isn't used much any more. The progression from this, when shades of grey became available, was to assign white to a negative amplitude (left of the zero line), through the shades of grey to mid grey at the zero line and on to black at a high positive amplitude (right of the zero line). The wave itself was no longer drawn, just the grey appropriate to the amplitude. This would give us the common method of B-Scan display used today, as shown in the example below.

B-Scan, current display method

Thus ends the history lesson. The example B-scan above, which I will use as an example throughout the rest of this blog post is from a town called Ewell. During the excavation by the Surrey Archaeological Society of a Roman road called Stane Street, I went to an adjacent modern road and walked along the pavement, crossing Stane Street. The curving feature in the middle is what I hope is the camber of the Roman road.

So what can you do with Snuffler? To get the data into Snuffler, there is a SEGY file import. SEGY is a common standard for sharing GPR data that can be used by many pieces of GPR software. Currently, Snuffler can currently only import a subset of this, the '2-byte, two's complement integer' type output by ReflexW. If you are trying this out and can't get Snuffler to import your file, zip it up and email it to me and I will get it working. The example above has already been processed in ReflexW, with static correction, background removal, gain and bandpass filters applied; basic stuff that Snuffler cannot do.

The standard method of display, as shown above, is the first thing you will see when you open the file in Snuffler, and here I want to discuss why I have a problem with this display method. There are two pieces of information that an image like this can tell you, amplitude and phase. It will show you both, but it will show you neither very well.

Amplitude is the strength of the wave. You can easily see a broad band across the top of the image with higher amplitude shown as lines of black and white rather than the grey of low amplitude. In a normal earth resistance plot in Snuffler, the scale might go from black as a low reading to white as a high reading. With GPR, both black and white are high amplitude, cutting the palette in half as far as contrast is concerned. It also doesn't help that the shades representing the high points on the wave only appear briefly before plunging back towards the zero line, making comparison of amplitude even more difficult. You can't even concentrate the display range at a certain amplitude to increase contrast without losing the other (negative) half of your data. The solution to this is quite common, but is not often applied to B-Scan data, the envelope. The envelope basically takes absolute value of the amplitude (makes everything positive) and fills in the gaps between the waves. This is normally done using a Hilbert Transform, but since my maths isn't that great, I literally took the absolute values and filled in the gaps instead. The result looks something like this. This, and the other display methods shown hereare not filters in the sense that they change the data, they are on the fly display methods that leave the data alone.

"Envelope" data

Now black is low amplitude and white is high amplitude. We have completely lost the phase information (more on that in a bit), but we can see the amplitude a lot better. We can still improve on this.
 Graph of amplitude data

The graph shows the amplitude on the X axis and the number of readings with that amplitude on the Y axis. If we move the low bound very slightly so the majority of the low readings are ignored and change the palette usage to non-linear, we can get better contrast. A non-linear palette is where instead of a shade of grey being assigned to an amplitude range, it is assigned to an equal number of readings with a similar amplitude. The result looks like this.

"Envelope" data with better contrast

So that's amplitude, what about phase? In simplified terms, a wave with a single frequency, given no obstacles, will have the same wavelength all along its length. There would be the same number of samples between the peaks of the wave. When the wave hits an obstacle, or in the case or GPR, hits a material with a different dielectric potential, the wavelength is temporarily shortened at the point of the interface before continuing, thus the phase of the wave changes. This effect is what produces the nice curve of the Roman road camber in our example. So what is the problem with the standard display method here? While it is easy to see the phase when the amplitude is high, the bottom of the image is a sea of grey with no contrast, rendering the phase there all but invisible. The common solution to this, often used in the oil industry where they go to very great depths, is something called Automatic Gain Control. For each line at a depth along the B-Scan, an average is taken in a time window around that line and the values along the line are changed to reflect the average. The effect of this is to remove the high amplitude bias at the top of the image and flatten it out, losing the amplitude information. This produced something like this.

Automatic Gain Control

Now you can see the lines all the way to the bottom, but we can still improve on this, as there isn't much contrast.


Graph of AGC data


If we change the display bounds to around the curve of the readings and use a non-linear display process as before, we can get a lot more contrast, resulting in a much clearer image as below.

Automatic Gain Control with better contrast

So far, so standard. Here is something that is hopefully a little different, but someone else has probably thought of it already. With the "envelope" data above, the high amplitude features are visible near the surface, but what about features at depth that have a high amplitude compared to other areas at the same depth, but a low amplitude compared to the features at the surface. They are invisible. If we apply automatic gain control followed by an envelope however, we get something like this.

AGC plus envelope

The effect is to lose the absolute amplitude we had before and gain a relative amplitude showing the strongest features at each sample depth. Of course, we can adjust the display parameters to get a better contrast.

Graph of AGC plus envelope data

If we set the low display bound to exclude the low amplitude samples and use a non-linear palette, we get this.

AGC plus envelope with better contrast

Now we can see the interesting features a lot clearer, including that interesting block to the left of Stane Street and Stane Street itself now shows as much more distinct than the rest of the material at that level. So we have three different display methods that are all better than the standard at doing their thing, but you need to look at all three to get the best picture, which is why most people stick with the standard display process that is a Jack of all trades but a master of none. Wouldn't it be great if you could see all three at once?

Channel merged image

Funky. Could you imagine scrolling through timeslices of that? It would be mindbending. Just like the channel merging I did a few versions back, this display method lets you view all three display methods in the same place. The "envelope" data is in the red channel, the automatic gain control in green and AGC plus envelope in blue. You lose some contrast having them all merged together, but the resulting image I think is far more useful than the standard display method that we are all used to. So there you go, itch scratched and funky GPR plots accomplished. I don't think I'll ever make Snuffler a fully functional GPR processing system, but hopefully this will be useful to someone out there. I may do time slices.

You can download this version of Snuffler in the usual place.

03 January 2016

It's talking Season Again

Winter is here, or not with these warm temperatures, so it is time for me to put down the geophys gear and do some talks.

The first will be for the Eastbourne Natural History and Archaeological Society. It is at St Saviour’s Church Hall, Spencer Rd, Eastbourne on Friday 8th of January 2016 at 7:30pm. I'll be doing some geophysics for the society in 2016, so I've been asked to talk about the sort of things that geophysics can find by talking about the settlements along the Hardham to Pevensey Roman road, focussing on Barcombe.

The second is at this year's symposium, run by the Sussex School of Archaeology, on the subject of the radar I did on Roman Chichester. I'll hopefully be back in 2016 to do more of that, so there may well be further talks. The date for that talk is the 19th of March, 2016, at the Huxley Building, University of Brighton.

08 November 2015

Processing UK Environment Agency LIDAR Data Tutorial

Recently, the Environment Agency has released its LIDAR data to the public. This is not the first time that DEM (Digital Elevation Map) data has been released from free. Satellite based DEM data has been released by ASTER and SRTM (SRTM is less noisy), but those have a horizontal resolution of 30m, which is quite coarse. The EA data, which was collected using LIDAR, has a much higher horizontal resolution, but the downside is that there is gaps in the data. Several of my archaeology friends have asked me to do tutorial on how to display this data.

If you look at the website, zoom in and click on a map tile, you can see that you can download 6 different sets of data. There is DTM (Digital Terrain Map) and DSM (Digital Surface Map) data. The difference is that the DTM data has stuff like trees stripped, which makes it much more useful for archaeology than the DSM data. It also comes in 3 resolutions, 0.5m, 1m & 2m. The trade off is that there is less area convered for the higher resolutions. I would recommend starting with the 2m data and downloading 1m for the same area if you really need it.

When you download a file, you will get a zip file containing up to 100 ASCII grid files, each covering a 1km square. This is not a normal image, but a file containing the height data. which we must process using GIS software. This tutorial will cover processing that data using a popular open source GIS package called QGIS, which you can download for free as a standalone application or as part of the OSGeo4w package. Here are the steps to follow.

If you haven't already got a project in QGIS, create a new one, setting the map projection to OSGB36 (EPSG:7405).

Unzip the contents of each zip file, and place the contents of each one in a folder of their own, for example, if your zip file is called LIDAR-DTM-2M-TR14.zip, create a directory called LIDAR-DTM-2M-TR14 and copy the contents of the zip file in there.

You will notice there are a lot of files, and we don't really want to deal with them all individually, so we will merge them together into one giant ASCII grid, covering a 10x10km square. To do this in QGIS, go to the Raster > Miscellaneous > Merge... menu. Next to the Input Files, press the Select... button and select all of the asc files in a directory. You also need an Output File for the merged data, for which I use the name of the directory again, which must be followed by a file extension of .asc, for example LIDAR-DTM-2M-TR14/LIDAR-DTM-2M-TR14.asc, preceeded of course by the rest of the file path to that directory. Finally, tick No data value and set it to -9999. Press ok and it will start working, which may take a while. It should load into QGIS automatically, having asked you what projection to use (OSGB36 again), but sometimes the process crashed. Don't worry, because the creation of the file has completed, so back on the QGIS main screen, press the Add Raster button and load the newly created file. Having merged all the files, you can delete all the original small unmerged files if you want, since they are pretty big. By default, you get the data displayed as a grey scale image with white as high and dark as low, as shown below.



You will notice that if you have more than one of these next to eachother, the edges don't match, so we need to set them to have a matching palette. We will also use a multi colour palette to provide a bit more contrast. Pick an image, right click it in the project window and select Properties. Go to the Style tab and change the Render type to Singleband pseudocolor. Then use the + sign to add bands. I currently use 0m-Black, 20m-Green, 40m-Brown, 60m-Orange, 80m-yellow and 100m-white. You can of course make up your own palette to suit your tastes. It should look something like this.


You can then right click on your image in the project window and select Styles > Copy Style and then right click on the rest of the images to do Styles > Paste Style to avoid entering that palette information all over again. This should give you something like this.


This shows the height information very well, but if you want to see small changes that might be created by archaeology, then we need to add a hillshade to this. To do this, go to Raster > Analysis > DEM. In the Input file, select one of your merged files. For the output file, use the same filename, except add -HS just before the .asc. You will notice that the resulting file will completely obscure the height data below. To fix this, right click on your newly created hillshade image, go to Properties, go to the Transparency tab and move the slider until it is at 50%. You should now be left with something like this. If you zoom in with QGIS, you will see a lot of the smaller features. Happy hunting!




29 September 2015

Green Waste, A Growing Problem

A relatively new, but increasing, bane for geophysicists and metal detectorists alike, is the practice of spreading green waste on fields. Part of the drive for increased recycling in society, which is good, this material comes from our gardens, but is rarely pure. Many households are not too fussed with what they will put in their bins, so plastic and metal will end up in the green waste. Farmers will buy this from the council and then use it to fertilise their fields. The metal component of this will cause problems for magnetometers, producing white noise from thousands of tiny dipoles. For example, this is a Roman settlement :


This is a Roman villa :


Once it's in there, it's there for good. Metal detectorists get it worse, as they are also picking up non-ferrous material. There is even a blog dedicated to the problem.