Processing Landsat 8 bands in Grass GIS

QGIS seems to get all the attention these days and perhaps rightly so. Grass on the other hand shows no sign of waning in the 'SAAS' world, if anything it has a continuing loyal bunch of users but in South Wales (UK) it is the outsider, it has an odd non microsoft like non-conforming interface. If you listen to the I.T departments it is too difficult for users to grasp, has an old fashioned command line and odd sounding module names like "i.fusion.brovey". Two years ago i would have walked the other way but since version 6 and a Python GUI there really is worthy alternative.

As with anything considered 'too' difficult for newcomers (fed up with 'newby') the long nights spent trying to grapple with 'mapsets' and 'regions' will eventually give way and reward you with the power and control over your work flows, customisations, alternative's and options in fact the Grass complexity hides a very simple no nonsense 'frill free' GIS. You could of course use a large number of the modules via a Web Feature Service in your own 'cloud' or other 'Geo' solution.

Landsat 8 (Operational Land Imager) launched on February 11th 2013 and images have been collected for some months now. These images are free and cover the earth on a sixteen day cycle there are eleven bands consisting of eight spectral bands 1 - 7, 'Band 8' is a panchromatic band (15m resolution) there is a new 'Band 1' useful for coastal/aerosol studies with 'Band 9' useful for Cirrus Cloud detection. Thermal Bands 10 and 11 measure surface temperature.

So enough talking and let's get started on this tutorial.

There are essentially only seven parts to this process and they are:

1.  Download tiff files from Earth Explorer
2.  Import the eleven bands into Grass GIS
3.  Initiate module 'i.landsat.toar'
4.  Initiate module 'r.colors' to 'histogram equalise' each band.
5.  Initiate module 'i.fusion.brovey' to pan sharpen bands 1 - 7
6.  Initiate module 'r.colors' to 'histogram equalise the resulting three  'brovey'  RGB bands.
7.  initiate module 'r.composite' to flatten the RGB bands into a single image

I want some regular images downloaded from the 'EarthExplorer' web site when ever there is a cloud free example. I will eventually want to create a python script to automate the process using the built in Grass Modeller but for now i will do it manually.

Once you login (register if not a member) you will be faced with this screen where you need to enter your address or use the 'WRS2' 'Path' and 'Row' method. For South Wales in the UK it is Path '204' and Row '24' once you've done that click on the 'Data Sets' button.

Next we need to access the 'Landsat Archive' listing and tick 'L8 OLI/TIRS' then click on 'Results'.

Now we have the most recent image on the top followed by the preceeding sixteen days image and so on.

On the icon bar next to each image click on the 'download' icon

Check the last option 'Level 1 GeoTIFF Data Product' then 'Select Download Option'

Click 'Download'

Save your 'gz' or 'tar.gz' somewhere you can access it. Download done!

Open Grass GIS (I'm using Version 6.4.3) then go to the 'File' menu and click on 'import raster data' then 'common formats import [r.in.gdal]' 

 

Unlike my screen shot above yours will contain a list of all the 11 bands ticked and ready for import under 'list of GDAL layers'.Once imported the layers will be listed in the Grass 'Layer Manager'.

 

If you tick 'LC82040242013363LGN00_B7@PERMANENT' you will notice that all Landsat bands are often very dark and raw prior to any adjustments making them about as useful as a chocolate spade in a heatwave.

 OK now we need to initiate the 'i.landsat.toar' module. Simply type in the module name in the 'Command console' tab then press 'return' on the keyboard.

Here is what it looks like.

Under the 'Required' tab enter 'LC82040242013363LGN00_B' in 'Base name' and 'LC82040242013363LGN00_refl' in 'Prefix', then in the 'metadata' tab we need to hit 'Browse' and locate the 'LC82040242013363LGN00_MTL.txt' file. This file contains all the band metadata saving you the need to enter any further values in this tab.

Next tab 'Settings' just leave it as it is.

Next tab 'Optional' again just leave it as it is.

Now click 'Run' then observe the following 'Command output' after about 13 minutes processing time.

 From here on we will only use these suffixed 'refl1', 'refl2', 'refl3' etc. If we open any of these new bands they will still look very dark at this stage. in order to remedy this we can apply a 'histogram equalisation' algorithm in order to spread out the most frequent intensity values. If you look at the image below this is what a band would look like before equalisation.

 And here it is after equalisation.

so we need to do this to all bands 1-7 or if you only need a natural looking final image simply equalise bands refl4, refl3 and refl2 corresponding to the Red, Green and Blue bands. In order to do this go back to the Grass 'Layer Manager' and click on 'Command Console' and type 'r.colors' and press 'return' on the keyboard.

Enter your required band (in the 'required' tab) as above using 'refl1' in this case. You will want to repeat this process for all the bands you need.

In the 'Colors' tab tick the 'Histogram equalisation' option and change the colour table to 'grey'.

Under the 'optional' tab leave it as it is. Now click on the 'Run' button.

 

I've run the 'r.colors' module on all bands.

Next we will initiate the 'i.fusion.brovey' module. Back to the Command Console and enter the module name as below and press the 'return' key.

This is what you should see.

One thing to note with this module and Landsat 8 images is the order in which you add the bands. The actual band combination for the above is refl4,refl3 and refl2 or normally known as 'Natural Colour' but the 'i.fusion.brovey' module was not designed with Landsat 8 bands in mind and the labelling is confusing. In Grass 7 there is a new 'i.pansharpen' module which will replace this. In the meantime add the 'refl' bands in reverse order refl2,refl3,refl4 and refl8 (pan sharp band), give a name to your output RGB bands such as 'brovey'.

In the 'Sensor' tab tick the 'Landsat Sensor' option.

No changes needed in the 'Optional' tab.

Note above you will now have three new RGB bands 'brovey.red', 'brovey.green' and 'brovey.blue' and the next stage is to equalise these bands again because the 'i.fusion.brovey' module seems to change the image back to it's original colour table (most annoying) perhaps it is possible to skip the original 'histogram equalization' module until we have initiated the 'i.fusion.brovey' module.

Next it's back to the Command Console to enter 'r.colors' and press the return key. After that you need to specify 'brovey.blue' as input raster map.

then under the 'Colors' tab and same as before tick 'Histogram equalisation' and the 'Grey' colour table.

Moving onto the 'Optional' tab which is left unchanged.

Once you've run this module on the remaining 'brovey.green' and brovey.red' we can have a quick look at the resulting image by using 'd.rgb'. This is a quick way to display all three RGB bands together without committing to creating a final image which will take approximately 30 minutes only to find your not going to be able to use that particular band combination. So, add your 'brovey.red', 'brovey.green' and 'brovey.blue' into the RGB selection drop downs below then press 'OK'.

The final RGB display will look something like this. Observe the purple urban areas the healthy vegetation in light greens and Industrial units in white and the sand and sediment in the shallow coastal waters in pale blue.

If your happy with this image you can take the final step and make a complete single band image using the 'r.composite' module.

Add your 'brovey' images in their respective RGB drop downs and add the final image name in this case I've added the number of the band, a description and band combination number as the 'output raster map'.

 

Under the 'Levels' tab leave all the default values unchanged.

Under the 'Optional' tab (no change to default).

OK that is it. You can now add your final single image into the view. I hope you find this tutorial useful. If your having difficulty or you think something is incorrect please let me know using the comment box.