Create Colour Composites for ALOS PALSAR Tiles

A common way to visualise dual-polarisation SAR data is to display a colour composite of HH (Red), HV (Green) and HH/HV(Blue).

Python Functions

Such composites can be created in RSGISLib using the following steps:

1. Create the ratio image using the bandMath.

   import rsgislib
   from rsgislib import imagecalc
   
   bandDefns = [imagecalc.BandDefn('hh', 'in_hh_file.kea', 1),
                imagecalc.BandDefn('hv', 'in_hv_file.kea', 1)]
   imagecalc.bandMath('out_hhhv_ratio.kea', 'hh/hv', 'KEA', \
                     rsgislib.TYPE_32FLOAT, bandDefns) 
   

2. Create the stack using stackImageBands

   import rsgislib
   from rsgislib import imageutils
   
   bands_list = ['in_hh_file.kea', 'in_hv_file.kea', 'out_hhhv_ratio.kea']
   band_names = ['HH','HV', 'HH/HV']
   gdaltype = rsgislib.TYPE_32FLOAT
   imageutils.stackImageBands(bands_list, band_names, 'out_hhhvratio_stack.kea', None, 
                                       0, 'KEA', rsgislib.TYPE_32FLOAT) 
   

3. Stretch using stretchImage

   import rsgislib
   from rsgislib import imageutils
   
   imageutils.stretchImage('out_hhhvratio_stack.kea, 'out_hhhvratio_stack_scale.kea', \
          False, '', True, False, 'GTiff', 
          rsgislib.TYPE_8INT, imageutils.STRETCH_LINEARSTDDEV, 2)
   

   This will create a tiff with values between 0 – 255 (8 bit), suitable for display in GIS packages or graphics programs such as Photoshop / GIMP.

Script for ALOS PALSAR tiles

For the 1×1 degree 25 m resolution ALOS PALSAR tiles, which can be freely downloaded from http://www.eorc.jaxa.jp/ALOS/en/palsar_fnf/fnf_index.htm for 2007 – 2010. A script (create_palsar_tiles_rgb.py) is available from RSGIS Scripts on Bitbucket. The script will also untar files, and can be run on downloaded data using:

python create_palsar_tiles_rgb.py —unzip downloaded_palsar_files

The script will untar each tile into a separate directory, create temp files for the ratio image and stack and save a tiff in with the original date files in the following structure:

N63W165_10_MOS
     |-----N63W165_10_MOS.tar.gz - Original tar.gz
     |-----N63W165_10_MOS_composite.tif - Composite generated
     |-----N63W165_10_date       - Acquisition date of each pixel
     |-----N63W165_10_date.hdr  
     |-----N63W165_10_linci      - Local incidence angle for each pixel
     |-----N63W165_10_linci.hdr 
     |-----N63W165_10_mask       - Data / no-data mask
     |-----N63W165_10_mask.hdr 
     |-----N63W165_10_sl_HH      - HH data (DN)
     |-----N63W165_10_s1_HH.hdr 
     |-----N63W165_10_sl_HV      - HV data (DN)
     |-----N63W165_10_s1_HH.hdr 

The resulting tiff will be something like the example below (Alaska; 63N165W)

Data Copyright JAXA

To apply the same stretch to multiple files a text file can be used with the minimum and maximum values for each band. For example:

#stddev
#band,img_min,img_max,out_min,out_max
1,48.4992,9286,0,255
2,396,4823.89,0,255
3,0.0181456,0.449104,0,255

This stretch file is specified using:

python create_palsar_tiles_rgb.py --unzip --stretchtxt stats.txt downloaded_palsar_files

It is also possible to create a stack using the coefficient of variance as the third band using:

# Coefficient of variance HH
python create_palsar_tiles_rgb.py --unzip --CoVHH downloaded_palsar_files 
# Coefficient of variance HV
python create_palsar_tiles_rgb.py --unzip --CoVHV downloaded_palsar_files

Obtaining JAXA PALSAR Data and Forest / Non-Forest Maps

JAXA have recently released their global forest / non-forest map at 50 m resolution and the Advanced Land Orbiting Satellite (ALOS) Phased Array L-Band SAR (PALSAR) data from which they were derived. This is really exciting because SAR data provides a different view of the world than optical data, which we’re more used to viewing. A particularly interesting feature of L-band SAR for mapping vegetation is the ability to ‘see’ through clouds and the canopy layer of vegetation. A good introduction to SAR data, in the context of vegetation mapping, is provided in the following paper:

Rosenqvist, A., Finlayson, C. M., Lowry, J and Taylor, D., 2007. The potential of long- wavelength satellite-borne radar to support implementation of the Ramsar Wetlands Convention. Aquatic Conservation: Marine and Freshwater Ecosystems. 17:229–244.

Downloading data

You can download data from:

http://www.eorc.jaxa.jp/ALOS/en/palsar_fnf/fnf_index.htm

You need to sign up for an account but this is a quick and straightforward process.

You can download data in 1 x 1 degree tiles or batches of 5 x 5 degrees. Data are in ENVI format, and can be read with GDAL, or programs that use GDAL (e.g., QGIS). If you don’t already have a viewer, you can download TuiView to open them with. ArcMap can read them (as it uses GDAL) but it won’t recognise it if you go through the ‘Add Data’ dialogue. However, you can just drag the files (larger files without the ‘.hdr’ extension) from windows explorer to the ‘Table of Contents’.

Mosaic data

To mosaic all files in a 5 x 5 degree batch (or any number of files), you can use a combination of GNU Parallel to untar and gdalbuildvrt. Assuming we want to mosaic the HH- and HV-polarisation PALSAR data the following commands can be used:

# Untar file
tar -xf N60W105_07_MOS.tar.gz

# Change into directory
cd N60W105

# Untar all files, in parallel using GNU Parallel
ls *.gz | parallel tar xf

# Create a list of HH and HV files
ls *_HH > hhfilelist.txt
ls *_HV > hvfilelist.txt

# Build VRT
gdalbuildvrt -input_file_list hhfilelist.txt N60W105_HH.vrt
gdalbuildvrt -input_file_list hvfilelist.txt N60W105_HV.vrt

This will create virtual rasters, which are text files containing references to the data. You can convert to real rasters (KEA, geotiff etc.,) using gdal_translate:

gdal_translate -of GTiff 60W105_HH.vrt 60W105_HH.tif
gdal_translate -of GTiff 60W105_HV.vrt 60W105_HV.tif

Calibrate data

The data are supplied as digital numbers, to calibrate and convert to dB, the following equation is used [1]:

You can run this calibration in RSGISLib using the following Python script:

 # Import RSGISLib
 import rsgislib
 from rsgislib import out image

# Set input and output image
 inimage = 'N60W105_HH.vrt'
 outimage = 'N60W105_HH_db.kea'

# Run image maths
 imagecalc.imageMath(inimage, outimage, '10*log10(b1^2) - 83.0,
         'KEA', rsgislib.TYPE_32FLOAT)

# Calculate stats and pyramids (for fast display)
 imageutils.popImageStats(outimage,True,0.,True)

Alternatively you could grab the calPALSAR_dB.py script to perform the same calculation using RIOS (which will run under windows, see earlier post for more details)

SAR data takes a while to get your head into but once you do it provides a wealth of information.

Update – data are now available at 25 m resolution from the same place

[1] Shimada, M. and Ohtaki, T. 2010. Generating large-scale high-quality SAR mosaic datasets: Application to PALSAR data for global monitoring. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. 3(4):637–656.