Tag Archives: hyperspectral

Working with hyperspectral data in ENVI BIL format using Python

For working with ENVI files I normally use GDAL as code can then be applied to different formats. However, there are a couple of limitations with GDAL when working with hyperspectral data in ENVI format:

  1. GDAL doesn’t copy every item from the header file to a new header file if they don’t fit in with the GDAL data model. Examples are FWHM and comments. Sometimes extra attributes are copied to the aux.xml file GDAL creates, these files aren’t read by ENVI or other programs based on IDL (e.g., ATCOR).
  2. For data stored Band Interleaved by Line (BIL) rather than Band Sequential (BSQ) reading and writing a band at a time is inefficient as it is necessary to keep jumping around the file

To overcome these issues NERC-ARF-DAN use their own Python functions for reading / writing header files and loading BIL files a line at a time. These functions have been tidied up and released through the NERC-ARF Tools repository on GitHub (https://github.com/pmlrsg/arsf_tools). The functions depend on NumPy.

To install them it is recommended to use the following steps:

  1. Download miniconda from http://conda.pydata.org/miniconda.html#miniconda and follow the instructions to install.
  2. Open a ‘Command Prompt’ / ‘Terminal’ window and install numpy by typing:
    conda install numpy
    
  3. Download ‘arsf_tools’ from GitHub (https://github.com/pmlrsg/arsf_tools/archive/master.zip)
  4. Unzip and within a ‘Command Prompt’ or ‘Terminal’ window navigate to the the folder using (for example):
    cd Downloads\arsf_tools-master
    
  5. Install the tools and library by typing:
    python setup.py install
    

Note, if you are using Linux you can install the arsf_binary reader from https://github.com/arsf/arsf_binaryreader which is written in C++. The ‘arsf_envi_reader’ module will import this if available as it is faster than the standard NumPy BIL reader.

If you are a UK based researcher with access to the JASMIN system the library is already installed and can be loaded using:

module load contrib/arsf/arsf_tools

An simple example of reading each line of a file, adding 1 to every band and writing back out again is:

from arsf_envi_reader import numpy_bin_reader
from arsf_envi_reader import envi_header

# Open file for output
out_file = open("out_file.bil", "w")

# Open input file
in_data = numpy_bin_reader.BilReader("in_file.bil")

for line in in_data:
out_line = line + 1
out_line.tofile(out_file)

# Copy header
envi_header.write_envi_header("out_file.bil.hdr",
                              in_data.get_hdr_dict())

# Close files
out_file.close()
in_data = None

A more advanced example is applying a Savitzky-Golay filter to each pixel. As the filter requires every band for each pixel it is efficient to work with BIL files.

For developing algorithms using spatial data, in particular multiple files it is recommended to convert the files to another band-sequential format using the ‘gdal_translate’ so they can be read using programs which use GDAL, for example RIOS or RSGISLib.

To convert files from ENVI BIL to ENVI BSQ and copy all header attributes the ‘envi_header’ module can be used after converting the interleave with GDAL. For example:

import os
import subprocess
from arsf_envi_reader import envi_header

# Set input and output image (could get with argparse)
input_image = "input.bil"
output_image = "output_bsq.bsq"

# Get interleave from file extension
output_interleave = os.path.splitext(output_image)[-1].lstrip(".")

# 1. Convert interleave
print("Converting interleave to {}".format(output_interleave.upper()))
gdal_translate_cmd = ["gdal_translate",
                      "-of", "ENVI",
                      "-co", "INTERLEAVE={}".format(output_interleave.upper())]
gdal_translate_cmd.extend([input_image, output_image])
subprocess.call(gdal_translate_cmd)

# 2. Copy header (GDAL doesn't copy all items)
# Find header files
input_header = envi_header.find_hdr_file(input_image)
output_header = envi_header.find_hdr_file(output_image)

# Read input header to dictionary
input_header_dict = envi_header.read_hdr_file(input_header)

# Change interleave
input_header_dict["interleave"] = output_interleave.upper()

# Write header (replace one generated by GDAL)
envi_header.write_envi_header(output_header, input_header_dict)

Useful GDAL commands for working with hyperspectral data

The Geospatial Data Abstraction Library (GDAL) contains a number of utilities for working with remote sensing data. The following are some commands which are particularly useful when working with airborne hyperspectral data.

Create subset and change format

To create a band and spatial subset of a larger hyperspectral dataset the following command can be used:

gdal_translate -of GTiff -co COMPRESS=LZW \
 -projwin 481500 766000 482500 765000 \
 -b 190 -b 40 -b 15  \
 outputs/f281013b_utm_wgs84N50_mapped_subsection.bil \
 outputs/f281013b_utm_wgs84N50_mapped_subsection_3band.tif
  • -of: Output Format – we are using GeoTiff file. For all available options see http://www.gdal.org/formats_list.html
  • -co: Creation Option – specifies the creation options, these are specific to the format. For GeoTiff we use compression to reduce the file size. For ENVI we can use the creation options to specify the Interleave (BSQ or BIL)
  • -projwin: A rectangle delineating an area to be extracted, in the format min X min Y max X max Y
  • -b: A given band to extract
  • final two arguments: input and output filenames

Create mosaic of multiple flightlines

Often a mosaic of all individual flightlines is required. Due to the large sizes of the mapped files this can often be time consuming. One option is to first make a ‘Virtial Raster’ using the ‘gdalbuildvrt’ command:

gdalbuildvrt -srcnodata 0 outputs/f2810_mapped_mosaic.vrt \
 outputs/f281013b_utm_wgs84N50_mapped_3band.bil \
 outputs/f281023b_utm_wgs84N50_mapped_3band.bil
  • -srcnodata: No data value of input lines
  • first argument: output VRT file
  • final arguments: list of input files

Rather than creating a new file a small .vrt text file is created which contains links to the individual lines.

Many programs which use GDAL (e.g., TuiView) can open this file directly. To make opening the file easier we can do some optimisation. First pre-calculate statistics using:

gdalinfo -stats outputs/f2810_mapped_mosaic.vrt
  • -stats: Calculate statistics and save to .aux.xml file
  • final argument: file name

The gdalinfo command can also be used to print information about a raster. Note for the .vrt file it shows the individual files which are included.

After creating statistics generate overviews (lower resolution copies at different resolutions) using gdaladdo:

gdaladdo -r average outputs/f2810_mapped_mosaic.vrt 2 4 8 16
  • -r average: Use average for resampling overviews
  • first argument: file name
  • numbers: overview levels to build (2x, 4x, 8x and 16x)

We can also make a real raster using gdal_translate (as shown previously).

Extract pixel values for a given location

gdallocationinfo -valonly -geoloc \
 outputs/f281013b_utm_wgs84N50_mapped_subsection.bil \
 481801 765289
  • -valonly: Only print values to the screen
  • -geoloc: Interprerate x and y values as coordinates rather than pixel and line
  • first argument: file name
  • x y: Easting and Northing of pixel to extract statistics from

This will print the values to the screen. To save as a text file we can redirect the output using ‘>’

gdallocationinfo -valonly -geoloc \
 outputs/f281013b_utm_wgs84N50_mapped_subsection.bil \
 481801 765289 > outputs/extracted_values.txt

This text file can be opened in other programs for further analysis.

This post is based on the hyperspectral practical session from the NERC Airborne Survey and Research Facility workshop held at Plymouth Marine Laboratory, March 2016. For more information about the ARSF Data Analysis Node (ARSF-DAN) and details about future training opportunities follow @ARSF_DAN on twitter. Hyperspectral data acquired by ARSF are available to download from NEODC.

Attribute LAS points using hyperspectral data

A good way of visualising LiDAR point clouds is to attribute them with colour information.
LAS files of point type 2 and above have the ability to store RGB colour information. These can be viewed with programs such as plas.io and CloudCompare.

To add RGB colour information to a LAS file ARSF-DAN have made a Python script available from: https://github.com/pmlrsg/arsf_tools

The script uses a combination of laspy (to read and write a LAS file) and GDAL (to read the raster) and is capable of working with images with more than 3 bands such as multispectral or hyperspectral files without creating a three band copy of the image first.

For this post I’ll demonstrate running the script on the JASMIN/CMEMS system. In a nutshell it is a Linux computing facility with lots of fast storage and direct access to many datasets, including the ARSF archive of airborne LiDAR and Hyperspectral data (apply for access to ARSF data here). Researchers in the UK who are part of the NERC or National Centre for Earth Observation (NCEO) scientific communities can apply for access. If you don’t have access to JASMIN you can run on your own machine just ignore the JASMIN specific steps.

  1. Apply for access and logon to JASMIN
    If you are using JASMIN the first stage is to get an account and log on by following the steps in this guide: http://jasmin.ac.uk/workflow/
  2. Set up scripts
    Once you have logged onto one of the shared VMs (e.g., cems-sci2.cems.rl.ac.uk) you need to load laspy, which isn’t installed by default. You can do this using:

    module load contrib/arsf/laspy/1.3.0
    

    Then checkout the repository containing the ‘colour_las_file.py’ script using:

    git clone https://github.com/pmlrsg/arsf_tools.git
    

    Check the script runs OK by typing:

    python2.7 arsf_tools/colour_las_file.py -h
    

    Note, you need to specify Python 2.7 as the default system Python is 2.6 but most of the Python libraries (e.g., GDAL) are build against 2.7. If you are not running the scripts on JASMIN just using ‘python’ should be fine.

  3. Colour using ARSF hyperspectral dataFor this example data from Mont Blanc flown as part of a EUFAR campaign in 2015 on Julian day 175a (24/06/2015) will be used. This data can be downloaded from NEODC (direct link).
    LIDAR_PATH=/neodc/arsf/2015/EUFAR15_28/EUFAR15_28-2015_175a_Mont_Blanc/LiDAR/flightlines/las1.2
    HS_PATH=/neodc/arsf/2015/EUFAR15_28/EUFAR15_28-2015_175a_Mont_Blanc/hyperspectral/flightlines/mapped
    
    python2.7 arsf_tools/colour_las_file.py \
              --image /vsizip/$HS_PATH/f175a183b_mapped_utm32n_bil.zip/f175a183b_mapped_utm32n.bil \
              --red 38 \
              --green 25 \
              --blue 12 \
              --scale \
             $LIDAR_PATH/LDR-EUFAR15_28-2015-175a-18.LAS \
             LDR-EUFAR15_28-2015-175a-18_col.LAS
    

    This will attribute using bands 38, 25 and 12 (true colour) and will scale the pixel values from 16 bit to 8 bit (0 – 255) using a standard deviation stretch.To read the hyperspectral data without unzipping a GDAL virtual file system is used, as described in a previous post.

  4. LASzip files and Download
    To reduce the size of the files before downloading you can use LASzip. On JASMIN LAStools can be loaded using:

    module load contrib/arsf/lastools/20150925
    

    To compress the LAS file and drop points flagged as noise the following command can be used:

    las2las -drop_class 7 -i LDR-EUFAR15_28-2015-175a-18_col.LAS -o LDR-EUFAR15_28-2015-175a-18_col.laz
    

    A similar command can be used to decompress files.

    See the following guide on transferring data off JASMIN: http://jasmin.ac.uk/how-to-use-jasmin/data-transfer/

  5. View using plas.io
  6. Open http://plas.io/ in a modern browser and click ‘Browse’ to load in the coloured LAS file.
    plasio_screenshot
    The nice thing about using plas.io and JASMIN is no specialist software needs to be installed on your local machine, just an ssh client and a web browser.

  7. View using Cloud Compare Viewer
  8. If you want to view files locally you can use Cloud Compare, which is available for Windows, OS X and Linux. To open the coloured LAS/LAZ file in the cloud compare viewer just drag it into the viewer window.

    cc_screenshot

  9. Additional – Colour using an existing three band image
    If you already have a three band image, with pixel values 0 – 255 a simpler command can be used:

    colour_las_file.py --image in_raster_file_rgb.tif \
                  in_las_file.las out_las_file.las
    

Copy information between ENVI header files

The header files (.hdr) of ENVI files are text files which contain the size of the file (pixels and lines) geospatial information and other meta data such as wavelengths, data units etc.,

Not all programs are aware of all items in then header file and some may become lost during processing. ARSF-DAN have made a script available which can copy attributes from one header to another. The script can be downloaded from https://github.com/pmlrsg/arsf_tools and is run using:

copy_header_info.py f221013_mapped_osng.bil.hdr \
                      f221013_nadir_mapped_osng.hdr \
                      -k "map info" "projection info" "wavelength"

To show the available fields the following script (also from the arsf_tools repo) can be used:

get_info_from_header.py f221013_mapped_osng.bil.hdr

Viewing hyperspectral data in TuiView

As TuiView is based on GDAL, it has always been able to open the ENVI files commonly used for hyperspectral data. However, there have been some recent enhancements to TuiView which have improved handling of hyperspectral data.

TuiView can be installed for Linux, OS X and Windows through conda using:

conda install -c rios tuiview

The source is available to download from https://bitbucket.org/chchrsc/tuiview/.

As hyperspectral data contain many bands, specifying the ones to use for display when you load the file in makes things easier. You can do this from the command line using the ‘–bands’ argument. For example:

tuiview --rgb --stddev --bands 280,460,160 f159223b_mapped_osng.bil

For the Specim Fenix dataset used as an example the command will display a colour composite of NIR, SWIR and Red at approximately the centre wavelengths of the Landsat 8 bands with a standard deviation stretch applied to the pixel values. Opening the image may take a while as the statistics need to be calculated for the the three bands selected for display in order to apply a standard deviation stretch.

Once the file has opened you can use the ‘Query tool’ to display spectral profiles for each pixel.

tuiview_screengrab

If you get an error that TuiView is unable to open the file it could be that ‘data ignore value = 0’ is set in the header and the GDAL driver is unable to handle it. You can export the environmental variable:

export GDAL_PAM_ENABLED=ON

on Linux/OS X, or:

set GDAL_PAM_ENABLED=ON

on Windows. This will create a ‘*.aux.xml’ metadata file with additional features the driver doesn’t support when TuiView opens the file. More information is available in the GDALPamDataset documentation.

If you have hyperspectral data which isn’t mapped (e.g., level1b / level2 data in the NASA Data processing levels) TuiView will not open these by default and will give an error about only allowing north-up images. If you know the data don’t have geospatial information you can force TuiView to open them by setting:

export TUIVIEW_ALLOW_NOGEO=YES

on Linux/OS X, or:

set TUIVIEW_ALLOW_NOGEO=YES

on Windows.

More details about available environmental variables used by TuiView are available on the TuiView Wiki.

Hyperspectral data from the NERC ARSF used for this post were provided courtesy of NERC via the NERC Earth Observation Data Centre (NEODC).