Robin's Blog

Behind the paper: Spatial variability of the atmosphere over southern England, and its effects on scene-based atmospheric corrections

I’ve just had my second journal paper published, and so I thought I’d start a series on my blog where I explain some of the background behind my publications, explain the implications/applications that my work has, and also provide a brief layman’s summary for non-experts who may be interested in my work. Hopefully this will a long-running series, with at least one post for each of my published papers – if I forget to do this in the future then please remind me!

So, this first post is about:

Wilson, R. T., Milton, E. J., & Nield, J. M. (2014). Spatial variability of the atmosphere over southern England, and its effect on scene-based atmospheric corrections. International Journal of Remote Sensing, 35(13), 5198-5218.

Publisher’s Link
Post-print PDF

Contributing areas for the Chilbolton AERONET site, calculated as part of the time-for-space substitution used to assess spatial variability

Layman’s summary

Satellite images are affected by the atmospheric conditions at the time the image was taken. These atmospheric effects need to be removed from satellite images through a process known as ‘atmospheric correction’. Many atmospheric correction methods assume that the atmospheric conditions are the same across the image, and thus correct the whole image in the same way. This paper investigates how much atmospheric conditions do actually vary across southern England, and tries to understand the effects of ignoring this and performing one of these uniform (or ‘scene-based’) atmospheric corrections. The results show that the key parameter is the Aerosol Optical Thickness (AOT) – a measure of the haziness of the atmosphere caused by particles floating in the air – and that it varies a lot over relatively small distances, even under clear skies. Ignoring the variation in this can lead to significant errors in the resulting satellite image data, which can then be carried through to produce errors in other products produced from the satellite images (such as maps of plant health, land cover and so on). The paper ends with a recommendation that, where possible, spatially-variable atmospheric correction should always be used, and that research effort should be devoted to developing new methods to produce high-resolution AOT datasets, which can then be used to perform these corrections.

Key conclusions

I always like my papers to answer questions in a way that actually affects what people do, and in this case there are a few key ‘actionable’ conclusions:

  1. Wherever possible, use a spatially-variable (per-pixel) atmospheric correction – particularly if your image covers a large area.
  2. Effort should be put into developing methods to retrieve high-resolution AOT from satellite images, as this data is needed to allow per-pixel corrections to be carried out.
  3. Relatively low errors in AOT can cause significant errors in atmospheric correction, and thus errors in resulting products such as NDVI. These errors may result from carrying out a uniform atmospheric correction when the atmosphere was spatially-variable, but they could just be due to errors in the AOT measurements themselves. Many people still seem to think that the NDVI isn’t affected by the atmosphere, but that is wrong: you must perform atmospheric correction before calculating NDVI, and errors in atmospheric correction can cause significant errors in NDVI.

Key results

  1. The range of AOT over southern England was around 0.1-0.5 on both days
  2. The range of PWC over southern England was around 1.5-3.0cm and 2.0-3.5cm on the 16th and 17th June respectively
  3. An AOT error of +/- 0.1 can cause a 3% error in the NDVI value

History & Comments

When I started my PhD I tried to find a paper like this one – and I couldn’t find one. I could find all sorts of comments in the literature – and in informal conversations with academics – that said that per-pixel atmospheric corrections were far better than scene-based corrections, but no-one seemed to have actually investigated the errors involved. So, I decided to investigate this myself as a sort of ‘Pilot Study’ for my PhD. This paper is basically a re-working of this Pilot Study.

Once I got started on this work I realised why no-one had done it! The first thing I needed to do was to use data on Aerosol Optical Thickness (AOT) and Precipitable Water Vapour (PWC) to find out how much spatial variation there is in these parameters. Unfortunately, the data is generally very low resolution, and so it is difficult to get a fair sense of how these parameters vary. In fact, almost half of the paper is taken up with describing the datasets that I’ve used, doing some validation on them, and then explaining how I used these datasets to estimate the range of values found over southern England during the dates in question. The datasets didn’t always agree particularly well, but we managed to establish approximate ranges of the values over the days in question.

Both of the days in question were relatively clear summer days, and I was surprised about the range of AOT and PWC values that we found. They were definitely nothing like uniform!

Once we’d established the range of AOT and PWC values, we performed simulations to establish the difference between a uniform atmospheric correction and a spatially-variable atmospheric correction. These simulations were carried out using Py6S: my Python interface to the 6S radiative transfer model. This made it very easy to perform multiple simulations at a range of wavelengths and with varying AOT and PWC values, and then process the data to produce useful results.

When performing a uniform atmospheric correction, a single AOT (or PWC) value is used across the whole image. We took this value to be the mean of the AOT (or PWC) values measured across the area, and then examined the errors that would result from correcting a pixel with this mean AOT when it actually had a higher or lower AOT. We performed simulations taking this higher AOT to be the 95th percentile of the AOT distribution, and the lower AOT to be the 5th percentile of this distribution. This meant that the errors found from the simulations would be found in at least 10% of the pixels in an image covering the study area.

Data, Code & Methods

Unfortunately I did the practical work for this paper before I had really taken on board the idea of ‘reproducible research’, so the paper isn’t easy to reproduce automatically. However, I do have the (rather untidy) code that was used to produce the results of the paper – please contact me if you would like a copy of this for any reason. The data are available from the following links – some of it is freely available, some only for registered academics:

  • MODIS: Two MODIS products were used, the MOD04 10km aerosol product and the MOD05 1km water vapour product. These were both acquired for tile h17v03 on the 16th and 17th June 2006, and are available to download through LADSWEB.
  • AERONET: Measurements from the Chilbolton site were used – available here.
  • GlobAerosol: Available via FTP from here.
  • Met Office visibility: Available to UK academics through BADC.
  • BIGF: Available to UK academics from BIGF.

What’s new in Py6S

The last few months have seen a flurry of activity in Py6S – probably caused by procrastinating from working on my PhD thesis! Anyway, I thought it was about time that I summarised the various updates and new features which have been released, and gave a few more details on how to use them.

These have all been released since January 2014, and so if you’re using version 1.3 or earlier then it’s definitely time to upgrade! The easiest way to upgrade is to simply run


pip install -U Py6S

in a terminal, which should download the latest version and get it all set up properly. So, on with the new features.

A wide range of bugfixes

I try to fix any actual bugs that are found within Py6S as soon as they are reported to me. The bugs fixed since v1.3 include:

  • More accurate specificiation of geometries (all angles were originally specified as integers, now they are specified as floating point values)
  • Fixed errors when setting custom altitudes in certain situations – for example, when altitudes have been set and then re-set
  • Fixes for ambiguity in dates when importing AERONET data – previously if you specified a date such as 01/05/2014 which could be interpreted either day-first (1st May) or month-first (5th January) then it assumed month-first, which was the opposite what the documentation specified. This now assumes day first – consistent with the documentation
  • Error handling has been improved for situations when 6S itself crashes with an error – rather than Py6S crashing it now states that 6S itself has encountered an error
  • Added the extraction of two outputs from the 6S output file that weren’t extracted previously: the integrated filter function and the integrated solar spectrum

Parallel processing support

Now when you use functions that run 6S for multiple wavelengths or multiple angles (such as the run_landsat_etm or run_vnir functions) they will automatically run in parallel. From the user’s point of view, everything should work in exactly the same way, but it’ll just be faster! How much faster, depends on your computer. If you’ve got a dual-core processor then it should be almost (but not quite) twice as fast. For a quad-core then it will probably be around three times faster, for an eight-core machine then it will probably be more like five times as fast. If you want to experiment then there is an extra parameter that you can pass to any of these functions to specify how many 6S runs to perform in parallel – just run something like:

run_landsat_etm(s, 'apparent_radiance', n=3)

to run three 6S simulations in parallel.

I’ve tested the parallel processing functionality extensively, and I’m very confident that it produces exactly the same answers as the non-parallel version. However, if you do run into any problems then please let me know immediately, and I’ll do whatever fixes are necessary.

Python 3 compatibility

Py6S is now fully compatible with Python 3. This has involved a number of changes to the Py6S source code, as well as doing some alterations to some of the dependencies so that they all work on Python 3 too. I don’t use Python 3 much myself, but all of the automated tests for Py6S now run on both Python 2.7 and Python 3.3 – so that should pick up any problems. However, if you do run into any issues, then please contact me.

Added wavelengths for two more sensors

Spectral response functions for Landsat 8 OLI and RapidEye are now included in the PredefinedWavelengths class, making it easy to simulate using these bands by code as simple as:

s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)

I’m happy to add the spectral response functions for other sensors – please email me if you’d like another sensor, and provide a link to the spectral response functions, and I’ll do the rest.

The future…

I’ve got lots of plans for the future of Py6S. Currently I’m finishing off my PhD, which is having to take priority over Py6S, but as soon as I’ve finished I should be able to release a number of new features.

Currently I’m thinking about ways to incorporate the building of Lookup Tables into Py6S – this should make running multiple simulations far quicker, and is essential to use Py6S for performing atmospheric corrections on images. I’m also considering a possible restructuring of the Py6S interface (or possibly a separate ‘modern’ Pythonic interface) for version 2.0 or 3.0. I’m also planning to apply to the Software Sustainability Institute Open Call, next year, with the aim of developing the software, and the community, further.


How to: Fix weird ENVI startup file issues

This post is more a note to myself than anything else – but it might prove useful for someone sometime.

In the dim and distant mists of time, I set up a startup file for ENVI which automatically loaded a specific image every time you opened ENVI. I have no idea why I did that – but it seemed like a good idea at the time. When tidying up my hard drive, I removed that particular file – and ever since then I’ve got a message each time I load ENVI telling me that it couldn’t find the file.

I looked in the ENVI preferences window, and there was nothing listed in the Startup File box (see below) – but somehow a file was still being loaded at startup. Strange.

ENVIConfig

 

I couldn’t find anything in the documentation about where else a startup file could be configured, and I searched all of the configuration files in the ENVI program folder just in case there was some sort of command in one of them – and I couldn’t find it anywhere.

Anyway, to cut a long story short, it seems that ENVI will automatically run a startup file called envi.ini located in your home directory (C:\Users\username on Windows, \home\username on Linux/OS X). This file existed on my machine, and contained the contents below – and deleting it stopped ENVI trying to open this non-existent file.

; envi startup script
open file = C:\Data\_Datastore\SPOT\SPOT_ROI.bsq

 


Simple parameter files for Python class-based algorithms

As part of my PhD I’ve developed a number of algorithms which are implemented as a class in Python code. An example would be something like this:

class Algorithm:
	def __init__(self, input_filename, output_basename, thresh, n_iter=10):
		self.input_filename = input_filename
		self.output_basename = output_basename
		
		self.thresh = thresh
		
		self.n_iter = n_iter
		
	def run(self):
		self.preprocess()
		
		self.do_iterations()
		
		self.postprocess()
		
	def preprocess(self):
		# Do something, using the self.xxx parameters

	def do_iterations(self):
		# Do something, using the self.xxx parameters
	
	def postprocess(self):
		# Do something, using the self.xxx parameters

The way you’d use this algorithm normally would be to instantiate the class with the required parameters, and then call the run method:

alg = Algorithm("test.txt", 0.67, 20)
alg.run()

That’s fine for using interactively from a Python console, or for writing nice scripts to automatically vary parameters (eg. trying for all thresholds from 0.1 to 1.0 in steps of 0.1), but sometimes it’d be nice to be able to run the algorithm from a file with the right parameters in it. This’d be particularly useful for users who aren’t so experienced with Python, but it can also help with reproducibility: having a parameter file stored in the same folder as your outputs, allowing you to easily rerun the processing.

For I while I’ve been trying to work out how to easily implement a way of using parameter files and the standard way of calling the class (as in the example above), without lots of repetition of code – and I think I’ve found a way to do it that works fairly well. I’ve added an extra function to the class which writes out a parameter file:

def write_params(self):
	with open(self.output_basename + "_params.txt", 'w') as f:
		for key, value in self.__dict__.iteritems():
			if key not in ['m', 'c', 'filenames']:
				if type(value) == int:
					valuestr = "%d" % value
				elif type(value) == float:
					valuestr = "%.2f" % value
				else:
					valuestr = "%s" % repr(value)

				f.write("%s = %s\n" % (key, valuestr))

This function is generic enough to be used with almost any class: it simply writes out the contents of all variables stored in the class. The only bit that’ll need modifying is the bit that excludes certain variables (in this case filenames, m and c, which are not parameters but internal attributes used in the class – in an updated version of this I’ll change these parameters to start with an _, and then they’ll be really easy to filter out).

The key thing is that – through the use of the repr() function – the parameter file is valid Python code, and if you run it then it will just set a load of variables corresponding to the parameters. In fact, the code to write out the parameters could be even simpler – just using repr() for every parameter, but to make the parameter file a bit nicer to look at, I decided to print out floats and ints separately with sensible formatting (two decimal places is the right accuracy for the parameters in the particular algorithm I was using – yours may differ). One of the other benefits of using configuration files that are valid Python code is that you can use any Python you want in there – string interpolation or even loops – plus you can put in comments. The disadvantage is that it’s not a particularly secure way of dealing with parameter files, but for scientific algorithms this isn’t normally a major problem.

The result of writing the parameter file as valid Python code is that it is very simple to read it in:

params = {}
execfile(filename, params)

This creates an empty dictionary, then executes the file and places all of the variables into a dictionary, giving us exactly what we’d want: a dictionary of all of our parameters. Because they’re written out from the class instance itself, any issues with default values will already have been dealt with, and the values written out will be the exact values used. Now we’ve got this dictionary, we can simply use ** to expand it to parameters for the __init__ function, and we’ve got a function that will read parameter files and create the object for us:

@classmethod
def fromparams(cls, filename):
	params = {}
	execfile(filename, params)
	del params['__builtins__']
	return cls(**params)

So, if we put all of this together we get code which automatically writes out a parameter file when a class is instantiated, and a class method that can instantiate a class from a parameter file. Here’s the final code, followed by an example of usage:

class Algorithm:
    def __init__(self, input_filename, output_basename, thresh, n_iter=10):
        self.input_filename = input_filename
        self.output_basename = output_basename
        
        self.thresh = thresh
        
        self.n_iter = n_iter

        self.write_params()

    def write_params(self):
        with open(self.output_basename + "_params.txt", 'w') as f:
            for key, value in self.__dict__.iteritems():
                if key not in ['m', 'c', 'filenames']:
                    if type(value) == int:
                        valuestr = "%d" % value
                    elif type(value) == float:
                        valuestr = "%.2f" % value
                    else:
                        valuestr = "%s" % repr(value)

                    f.write("%s = %s\n" % (key, valuestr))
            
    def run(self):
        self.preprocess()
        
        self.do_iterations()
        
        self.postprocess()

    @classmethod
    def fromparams(cls, filename):
        params = {}
        execfile(filename, params)
        del params['__builtins__']
        return cls(**params)
        
    def preprocess(self):
        # Do something, using the self.xxx parameters

    def do_iterations(self):
        # Do something, using the self.xxx parameters
    
    def postprocess(self):
        # Do something, using the self.xxx parameters

And the usage goes something like:

# Create instance with code
alg = Algorithm("input.txt", "output", 0.25, n_iter=20)
alg.run()

# Create instance from parameter file
alg = Algorithm.fromparams('output_params.txt')

How to: Find out what modules a Python script requires

I do a lot of my academic programming in Python, and – even though I often write about the importance of reproducible research – I don’t always document my code very well. This sometimes leads to problems where I have some code running fine, but I don’t know which modules it requires. These could be external libraries, or modules I’ve written myself – and it’s very frustrating to have to work out the module requirements by trial and error if I transfer the code to a new machine.

However, today I’ve realised there’s a better way: the modulefinder module. I’ve written a short piece of code which will produce a list of all of the ‘base’ or ‘root’ modules (for example, if you run from LandsatUtils.metadata import parse_metadata, then this code will record LandsatUtils) that your code uses, so you know which you need to install.

Hopefully, like me, it’ll save you some time.


How to: Solve GDAL error ‘An error occurred while writing a dirty block’

When running GDAL on my university’s supercomputer yesterday I got the following error:

ERROR 1: Landsat_Soton.tif, band 1: An error occured while writing a dirty block

This post is really just to remind me how to solve this error – I imagine the error may have a multitude of possible causes. In my case though, I knew I’d seen it before – and fixed it – but I couldn’t remember how. It turns out that it’s really simple: GDAL is giving an error saying that it can’t write part of the output file to the hard drive. In this case, it’s because the supercomputer that I’m using has quotas for the amount of storage space each user can use – and I’d gone over the quota ‘hard limit’, and therefore the operating system was refusing to write any of my files.

So, the simple answer is to delete some files, and then everything will work properly!

(If you’re not using a shared computer with quotas, then this may be because your hard drive is actually full!)

 

 


How effective is my research programming workflow? The Philip Test – Part 4

10. Can you re-generate any intermediate data set from the original raw data by running a series of scripts?

It depends which of my projects you’re talking about. For some of my nicely self-contained projects then this is very easy – everything is encapsulated in a script or a series of scripts, and you can go from raw data, through all of the intermediate datasets, to the final results very easily. The methods by which this is done vary, and include a set of Python scripts, or the use of the ProjectTemplate package in R. Since learning more about reproducible research, I try to ‘build in’ reproducibility from the very beginning of my research projects. However, I’ve found this very difficult to add to a project retrospectively – if I start a project without considering this then I’m in trouble. Unfortunately, a good proportion of my Phd is in that category, so not everything in the PhD is reproducible. However, the main algorithm that I’m developing is – and that is fully source-controlled, relatively well documented and reproducible. Thank goodness!

11. Can you re-generate all of the figures and tables in your research paper by running a single command?

The answer here is basically the same as above: for some of my projects definitely yes, for others, definitely no. Again, there seems to be a pattern that smaller more self-contained projects are more reproducible – and not all figures/tables of my PhD thesis can be reproduced – but generally you’ve got a relatively good chance. At the moment I don’t use things like Makefiles, and don’t write documents with Sweave, KnitR or equivalents – so to reproduce a figure or table you’ll often have to find a specific Python file and run it (eg. create_boxplot.py, or plot_fig1.py), but it should still produce the right results.

12. If you got hit by a bus, can one of your lab-mates resume your research where you left off with less than a week of delay?

Not really not – it would be difficult, even for my supervisor or someone who knew a lot about what I was doing to take over my work. My “bus factor” is definitely 1 (although I hope that the bus factor for Py6S is fractionally greater than 1). Someone who had a good knowledge of Python programming, including numpy, scipy, pandas and GDAL, would have a good chance at taking over one of my better-documented and more-reproducible smaller projects – but I think someone would struggle to pick up my PhD. In many ways though, that’s kinda the point of a PhD – you’re meant to end up being the World Expert in your very specific area of research, which would make it very difficult for anyone to pick up anyone’s PhD project.

For one of my other projects, it may take a while to get familiar with it – but it should be perfectly possible to take my code, along with drafts of papers and/or other documentation I’ve written and continue the research. In many ways that is the whole point of reproducible research: aiming to develop research that someone else can easily reproduce and extend. The only difference is that usually the research is reproduced/extended after it’s been completed by you, whereas if you get hit by a bus then it’ll never have been completed in the first place!


How to: Create a Landsat Metadata database, so you can select images by various criteria

Recently I ran into a situation where I needed to select Landsat scenes by various criteria – for example, to find images over a certain location, within a certain date range, with other requirements on cloudiness and so on. Normally I’d do this sort of filtering using a tool like EarthExplorer, but I needed to do this for about 300 different sets of criteria – making an automated approach essential.

So, I found a way to get all of the Landsat metadata and import it into a database so that I could query it at will and get the scene IDs for all the images I’m interested in. This post shows how to go about doing this – partly as a reference for me in case I need to do it again, but hopefully other people will find it useful.

So, to start, you need to get the Landsat metadata from the USGS. On this page, you can download the metadata files for each of the Landsat satellites separately (with Landsat 7 metadata split into SLC-on and SLC-off).

1_Website

You’ll want the CSV files, so click the link and have a break while it downloads (the CSV files are many hundreds of megabytes!). If you look at the first line of the CSV file once you’ve downloaded it (you may not want to load it in a text editor as it is such a huge file, but something like the head command will work fine), you’ll see the huge number of column headers giving every piece of metadata you could want! Of course, most of the time you won’t want all of the metadata items, so we want to extract just the columns we want.

The problem with this is that lots of the traditional tools that are used for processing CSV files – including text editors, database import tools and Excel – really don’t cope well with large CSV files. These Landsat metadata files are many hundreds of megabytes in size, so we need to use a different approach. In this case, I found that the best approach was using one of the tools from csvkit, a set of command-line tools for processing CSV files, written in Python. One of the key benefits of these tools is that they process the file one line at a time, in a very memory-efficient way, so they can work on enormous files very easily. To extract columns from a CSV file we want to use csvcut, and we can call it with the following command line

csvcut -c 5,6,1,2,3,7,8,20,25,30 LANDSAT_ETM.csv > LANDSAT_ETM_Subset.csv

This will extract the 5th, 6th, 1st, 2nd 3rd etc columns from LANDSAT_ETM.csv to LANDSAT_ETM_Subset.csv. To get a list of the columns in the file along with their ID number, so that you can choose which ones you want to extract, you can run:

csvcut -n LANDSAT_ETM.csv

After doing this you’ll have a far smaller CSV file in LANDSAT_ETM_Subset.csv that just contains the columns you’re interested in. There’s only one problem with this file – it still has the headers at the beginning. This is great for a normal CSV file, but when we import it into the database we’ll find that the header line gets imported too – not what we want! The easiest way to remove it is using the following command:

cat LANDSAT_ETM_Subset.csv | sed "1 d" > LANDSAT_ETM_Subset.csv

Again, this doesn’t load the whole file in to memory, so will work with large files happily.

We then need to create the database. This can be done with any database system, but to get a simple local database I decided to use SQLite. Once you’ve installed this you can do everything you need from the command-line (you can create the tables using a GUI tool such as SQLite Administrator, but you won’t be able to do the import using that tool – it’ll crash on large CSV files). To create a database simply run:

sqlite LandsatMetadata.sqlite

which will create a database file with that name, and then drop you in to the SQLite console. From here you can type any SQL commands (including those to create or modify tables, plus queries), plus SQLite commands which generally start with a . In this case, we need to create a table for the various columns we’ve chosen from our CSV. It is important here to make sure that the column names are exactly the same as those in the CSV, or the import command won’t work (you can change the names later with ALTER TABLE if needed). You can take the following SQL and modify it to your needs.

CREATE TABLE [images] (
[browseAvailable] BOOLEAN  NULL,
[browseURL] VARCHAR(500)  NULL,
[sceneID] VARCHAR(100)  NULL,
[sensor] VARCHAR(20)  NULL,
[acquisitionDate] DATE  NULL,
[path] INTEGER  NULL,
[row] INTEGER  NULL,
[cloudCoverFull] FLOAT  NULL,
[dayOrNight] VARCHAR(10)  NULL,
[sceneStartTime] VARCHAR(30)  NULL
);

Just type this into the SQLite console and the table will be created. We now need to import the CSV file, and first we have to define what is used as the separator in the file. Obviously, for a CSV file, this is a comma, so we type:

.separator ,

And then to actually import the CSV file we simply type:

.import LANDSAT_ETM_Subset.csv images

That is, .import followed by the name of the CSV file and the name of the table to import into. Once this is finished – it may take a while – you can check that it imported all of the rows of the CSV file by running the following query to get the number of rows in the table:

SELECT COUNT() from images;

and you can compare that to the output of

wc -l LANDSAT_ETM_Subset.csv

which will count the lines in the original file.

Your data is now in the database and you’re almost done – there’s just one more thing to do. This involves changing how the dates and times are represented in the database, so you can query them easily. Still in the SQLite console, run:

UPDATE images
SET startTime=time(substr(images.sceneStartTime,10, length(images.sceneStartTime)));

And then…you’re all done! You can now select images using queries like:

SELECT * FROM images WHERE path=202 AND row=24
AND acquisitionDate > date("2002-03-17","-1 months")
AND acquisitionDate < date("2002-03-17","+1 months")

 Once you’ve got the results from a query you’re interested in, you can simply create a text file with the sceneIDs for those images and use the Landsat Bulk Download Tool to download the images.


How effective is my research programming workflow? The Philip Test – Part 3

7. Do you use version control for your scripts?

Yes, almost always. I’ve found this a lot easier since I started using Git – to start using version control with Git simply requires running “git init” – whereas with SVN you had to configure a new repository and do all sorts of admin work before you could start. Of course, if you want to host the Git repository remotely then you’ll need to do some admin, but all of that can be done at a later date, once you’ve started recording changes to your code. This really helps me start using version control from the beginning of the project – there’s less to interrupt me when I’m in the mood to write code, rather than to fiddle around with admin.

Occasionally I still forget to start using version control at the beginning of a project, and I always kick myself when I don’t. This tends to be when I’m playing around with something that doesn’t seem important, and won’t have any future…but it’s amazing how many of those things go on to become very important (and sometimes get released as open-source projects of their own). I’m trying to remember to run git init as soon as I create a new folder to write code in – regardless of what I’m writing or how important I think it’ll be.

8. If you show analysis results to a colleague and they offer a suggestion for improvement, can you adjust your script, re-run it, and produce updated results within an hour?

Sometimes, but often not – mainly because most of my analyses take more than an hour to run! However, I can often do the alterations to the script/data and start it running within the hour – and colleagues are often surprised at this. For example, someone I’m currently working with was very surprised that I could re-run my analyses, producing maximum values rather than mean values really easily – in fact, within about five minutes! Of course, all I needed to do in my code was change the mean() function to the max() function, and re-run: simple!

I’ve already talked a bit on this blog about ProjectTemplate – a library for R that I really like for doing reproducible research, and the benefits of ProjectTemplate are linked to this question. By having defined separate folders for library functions and analysis functions – and by caching the results of analyses to stop them having to be run again (unless they’ve changed), it really helps me know what to change, and to be able to produce new results quickly. I’d recommend all R programmers check this out – it’s really handy.

9. Do you use assert statements and test cases to sanity check the outputs of your analyses?

Rarely – except in libraries that I develop, such as Py6S. I recently wrote a post on how I’d set up continuous integration and automated testing for Py6S, and I am strongly aware of how important it is to ensure that Py6S is producing accurate results, and that changes I make to improve (or fix!) one part of the library don’t cause problems elsewhere. However, I don’t tend to do the same for non-library code. Why?

Well, I’m not really sure. I can definitely see the benefits, but I’m not sure whether the benefits outweight the time that I’d have to put in to creating the tests. Assert statements would be easier to add, but I’m not sure if I should start using them in my code. Currently I use exception handling in Python to catch errors, print error messages and – crucially – set all of the results affected by that error to something appropriate (normally NaN), and these have the benefit that the code continues to run. So, if one of the 100 images that I’m processing has something strange about it then my entire processing chain doesn’t stop, but the error gets recorded for that image, and the rest of the images still get processed. I actually think that’s better – but if someone has an argument as to why I should be using assert statements then I’d love to hear it.


Software choices in remote sensing

I recently read the article Don’t be a technical masochist on John D. Cook’s blog, and it struck a chord with me about the way that I see people choosing software and programming tools in my field.

John states “Sometimes tech choices are that easy: if something is too hard, stop doing it. A great deal of pain comes from using a tool outside its intended use, and often that’s avoidable…a lot of technical pain is self-imposed. If you keep breaking your leg somewhere, stop going there.”.

I like to think that I can do some GIS/remote-sensing analyses more efficiently than others – and I think this is because I have a broad range of skills in many tools. If the only GIS/Remote Sensing tool you know how to use is ArcGIS then you try and do everything in Arc – even if the task you’re doing is very difficult (or even impossible) to do within Arc. Similarly, if you only know Python and won’t touch R, then you can’t take advantage of some of the libraries that are only available in R, which might save you a huge amount of time. I wouldn’t say I’m an expert in all of these tools, and I prefer some to others, but I have a working knowledge of most of them – and am able to acquire a more in-depth knowledge when needed.

Don’t get me wrong, sometimes it is very important to be able to do things within a certain tool – and sometimes it’s worth pushing the boat out a bit to try and see whether it’s possible to get a tool to do something weird but useful. Often though, it’s better to use the best tool for the job. That’s why, if you watch me work, you’ll see me switching backwards and forwards between various tools and technologies to get my job done. For example:

  • I’m currently working on a project which involves a lot of time series manipulation. When I started this project, the Python pandas library (which deals very nicely with time series) wasn’t very mature, and I wasn’t prepared to ‘bet the project’ on this very immature library. So, even though Python is my preferred programming language, I chose to use R, with the xts library to handle my time-series analysis.
  • I don’t use ArcGIS for remote sensing analysis, and I don’t use ENVI to do GIS. Yes, both programs allow you to deal with raster and vector data, but it’s almost always easier to use ENVI to do the remote sensing work, and then transfer things into ArcGIS to overlay it with vector data or produce pretty output maps (ENVI’s map output options are pretty awful!). I’ve lost track of the number of students I’ve found who’ve been really struggling to do satellite data processing in ArcGIS that would have taken them two minutes in ENVI.
  • If there’s a great library for a different programming language then I use it. For example, I recently needed to create a set of images where each pixel contained a random value, but there was spatial correlation between all of the values. I investigated various tools which purported to be able to do this (including some random Python code I found, ArcGIS, specialist spatial stats tools and R) and in the end the one that I found easiest to get it to work in was R – so that’s what I used. Yes, it meant I couldn’t drive it easily from Python (although using RPy, or, as a last resort, running R from the command line using the Python subprocess module) but that was a far easier option compared to trying to write some code to do this from scratch in Python.

Overall, the tools I use for my day-to-day work of GIS/Remote-sensing data processing include (in approximate order of frequency of use): ENVI, GDAL’s command-line tools, QGIS, ArcGIS, GRASS, eCognition, PostGIS, Google Earth, BEAM, and probably more that I can’t think of at the moment. On top of that, in terms of programming languages, I use Python the most, but also R and IDL fairly frequently – and I’ve even been known to write Matlab, Mathematica and C++ code when it seems to be the best option (for example, Mathematica for symbolic algebra work).

Having a basic knowledge of all of these tools (and, of course, having them installed and set up on my computers and servers) allows me to get my work done significantly faster, by using the best (that is, normally, easiest) tool for the job.