Worked example updated

docs/source/worked_example.rst

@@ -3,11 +3,13 @@
 
 Here we'll run through an example based on an area of Southeastern Tanzania called Kilwa District (Fig). This is a location where the University of Edinburgh coordinates a network of forest plots where we have good data on aboveground biomass.
 
+We'll download data, generate maps of gamma0 backscatter, AGB and forest cover.
+
 
 Downloading data
 ----------------
 
-To follow.
+[To follow]
 
 Open Python and import biota
 ----------------------------
@@ -33,8 +35,8 @@ To load the biota module, type:
     
     >> import biota
 
-Loading a tile
---------------
+Loading an ALOS tile
+--------------------
 
 Data from the ALOS mosaic is provided as a series of 1 x 1 degree tiles. To load a tile in memory, we need to tell ``biota`` what directory the ALOS mosaic data are stored in, and what latitude and longitude we want to load. To save us from writing them out repeatedly, we can store these as variables:
 
@@ -69,7 +71,7 @@ The new object called ``tile_2007`` has a range of attributes. These can be acce
     >>> tile_2007.xRes, data_2007.yRes # Pixel resolution in meters
     (24.401, 24.579)
 
-*Advanced:* The tile also contains projection information for interaction with ``GDAL``:
+**Advanced:** The tile also contains projection information for interaction with ``GDAL``:
 
 .. code-block:: python
     
@@ -80,8 +82,8 @@ The new object called ``tile_2007`` has a range of attributes. These can be acce
     >>> tile_2007.proj # Projection wkt
     'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
 
-There are a range of other options that can be specified when opening an ALOS tile, but we'll return to these.
-    
+There are a few other options that can be specified when loading an ALOS tile, but we'll return to these in the see the :ref:`furtheroptions` section.
+
 Extracting backscatter information
 ----------------------------------
 
@@ -138,40 +140,163 @@ If we want to save this data to a geoTiff, we can run:
     
     >>> gamma0_2007 = tile_2007.getGamma0(polarisation = 'HV', units = 'decibels', output = True)
 
-which will write a GeoTiff file to the current working directory. This file can be processed and visualised in standard GIS and remote sensing software (e.g. QGIS, GDAL).
+This will write a GeoTiff file to the current working directory. This file can be processed and visualised in standard GIS and remote sensing software (e.g. QGIS, GDAL).
 
 Building a map of AGB
 ---------------------
 
-In a similar way to loading gamma0 backscatter, we can generate images of AGB. Note: by default ``biota`` uses an equation calibrated for ALOS-1 in miombo woodlands of Southern Africa. It's advisable to have a locally calibrated biomass-backscatter equation to improve on predictions.
+In a similar way to loading gamma0 backscatter, we can show maps of AGB.
 
 .. code-block:: python
 
-        agb_2007 = tile_2007.getAGB(show = True) # To display AGB
-        agb_2007 = tile_2007.getAGB(output = True) # To output AGB map to a GeoTiff
+    >>> agb_2007 = tile_2007.getAGB(show = True)
+
+Areas in darker green show denser forest:
 
 .. figure:: images/agb.png
 
+Like the previous function (and most others in the ``biota`` module), we can output a GeoTiff as follows:
+
+.. code-block:: python
+
+    >>> agb_2007 = tile_2007.getAGB(output = True) # To output AGB map to a GeoTiff
+
+.. note:: By default ``biota`` uses an equation calibrated for ALOS-1 in miombo woodlands of Southern Africa. It's advisable to have a locally calibrated biomass-backscatter equation to improve predictions.
+
 Building a forest cover map
 ---------------------------
 
-With a few modifications to the above scipt, we can generate a map of forest cover baed on thresholds of minimum AGB and minimum area. Here we assume a forest definition consistent with a minimum AGB of 10 tC/ha over a minimum area of 1 hectare.
+A forest cover map (or 'woody cover') can be generated as follows:
+
+.. code-block:: python
+
+    >>> woodycover_2007 = tile_2007.getWoodyCover(show = True)
+
+and output:
+
+.. code-block:: python
+
+    >>> woodycover_2007 = tile_2007.getWoodyCover(output = True)
+
+.. figure:: images/woodycover.png
+
+By default ``biota`` will use a generic definition of forest of 10 tC/ha with no minimum area. In the next section we'll discuss how this can be customised.
+
+Further options when loading an ALOS tile
+-----------------------------------------
+.. _furtheroptions:
+
+``biota`` supports a range of options for data processing and forest definitions. These options should be specified when loading a tile.
+
+Speckle filtering
+~~~~~~~~~~~~~~~~~
+
+Radar data are often very noisy as the result of 'radar speckle', which can be supressed with a speckle filter. The ``biota`` module has an Enhanced Lee speckle filter, which can be applied to the ALOS tile by loading it as follows:
+
+.. code-block:: python
+    
+    >>> tile_2007 = biota.LoadTile(data_dir, latitude, longitude, 2007, lee_filter = True)
+
+Downsampling
+~~~~~~~~~~~~
+
+Data volumes can be reduced through downsampling. This comes at a cost to resolution, but does have the positive effect of reducing speckle noise. For example, to halve the resolution of output images, set the parameter ``downsample_factor`` to 2:
+
+.. code-block:: python
+    
+    >>> tile_2007 = biota.LoadTile(data_dir, latitude, longitude, 2007, downsample_factor = 2)
+
+Changing forest definitions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+For many purposes it's useful to classify regions into forest and nonforest areas. To achieve this with ``biota`` a threshold AGB (``forest_threshold``) and a minimum area (``area_threshold``) that separate forest from nonforest can be specified. For example, for a forest definition of 15 tC/ha with a minimum area of 1 hecatare:
+
+.. code-block:: python
+    
+    >>> tile_2007 = biota.LoadTile(data_dir, latitude, longitude, 2007, forest_threshold = 15, area_threshold = 1)
+
+Changing output directory
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The current working directory may not be the best place to output GeoTiff files. An output directory can be specified as follows:
+
+.. code-block:: python
+    
+    >>> tile_2007 = biota.LoadTile(data_dir, latitude, longitude, 2007, output_dir = '~/outputs/)
+
+Masking data
+------------
+
+[To follow]
+
+Putting it all together
+-----------------------
+
+Using the commands above, we can create a script to automate the pre-processing of an ALOS tile to generate outputs of gamma0 (HV backscatter in units of decibels), AGB and forest cover for the year 2007. We'll filter the data and specify a forest threshold of 15 tC/ha with a minimum area of 1 hectare, Using a text editor:
 
 .. code-block:: python
+    
+    # Import the biota module
+    import biota
+
+    # Define a variable with the location of ALOS tiles
+    data_dir = '~/DATA/'
+    
+    # Define and output location
+    output_dir = '~/outputs/'
+    
+    # Define latitude and longitude
+    latitude = -9
+    longitude = 38
+    
+    # Load the ALOS tile with specified options
+    tile_2007 = biota.LoadTile(data_dir, longitude, latitude, 2007, lee_filter = True, forest_threshold = 15., area_threshold = 1, output_dir = output_dir)
+    
+    # Calculate gamma0 and output to GeoTiff
+    gamma0_2007 = tile_2007.getGamma0(output = True)
+    
+    # Calculate AGB and output to GeoTiff
+    gamma0_2007 = tile_2007.getAGB(output = True)
 
+    # Calculate Woody cover and output to GeoTiff
+    gamma0_2007 = tile_2007.getWoodyCover(output = True)
+    
+Save this file (e.g. ``process_2007.py``), and run on the command line:
+
+.. code-block::
+    
+    python process_2007.py
+
+**Advanced:** To process multiple tiles, we can use nested ``for`` loops. We can also add a ``try``/``except`` condition to prevent the program from crashing if an ALOS tile can't be loaded (e.g. over the ocean).
+
+.. code-block:: python
+    
     # Import the biota module
     import biota
 
     # Define a variable with the location of ALOS tiles
-    data_dir = '~/SMFM/ALOS_data/'
-
-    # Loop through the latitudes/longitudes that cover Kilwa District
-    for longitude in range(35, 40):
-        for latitude in range(-10, -5):
-                    
-            # Load the tile for the year 2010. Apply a speckle filter (lee_filter = True).
-            # Apply thresholds for forest cover (10 tC/ha) with a minimum area (1 hectare)
-            tile = biota.LoadTile(data_dir, longitude, latitude, 2010, lee_filter = True, forest_threshold = 10., area_threshold = 1)
+    data_dir = '~/DATA/'
+    
+    # Define and output location
+    output_dir = '~/outputs/'
+    
+    for latitude in range(35, 40):
+        for longitude in range(-10,-5):
+        
+            # Load the ALOS tile with specified options
+            try:
+                tile_2007 = biota.LoadTile(data_dir, longitude, latitude, 2007, lee_filter = True, forest_threshold = 15., area_threshold = 1, output_dir = output_dir)
+            
+            except:
+                continue
             
-            # Calculate forest cover and output to a GeoTiff
-            tile.getWoodyCover(output = True)
+            # Calculate gamma0 and output to GeoTiff
+            gamma0_2007 = tile_2007.getGamma0(output = True)
+            
+            # Calculate AGB and output to GeoTiff
+            gamma0_2007 = tile_2007.getAGB(output = True)
+
+            # Calculate Woody cover and output to GeoTiff
+            gamma0_2007 = tile_2007.getWoodyCover(output = True)
+    
+