Radiometric Normalization

This library implements functionality for normalizing a candidate image to time-invariant features in a set of reference images covering a time series. This includes generating a time-invariant reference image from a time stack, identifying features that are invariante between the reference time stack and a candidate image, calculating a linear transformation that normalizes the candidate image to the reference image, applying the linear transform, and validating the results.

The primary use case for this library is radiometrically normalizing a satellite image to a time series from a reference dataset.

Development Environment

This library uses a Vagrant VM for the development environment and requires Vagrant on the host computer.

To start the VM, in the root directory of the repo type:

vagrant up

Log into the VM by typing:

vagrant ssh

Navigate to the root directory:

cd /vagrant

This directory is shared between the VM and host.

Run the unit tests by typing:

nosetests ./tests

Organization of the repo

The code in this repo is kept in two directories: 'radiometric_normalization' and 'tests'

'radiometric_normalization' contains the algorithm and functions of the library:

'time_stack.py' is the module that can calculate a time stack (average image) from a series of images (depreciated).
'pif.py' contains the functions that find pseudo invariant feature pixels on within an image.
'transformation.py' calculates the transformations to make the candidate data more consistent with the reference data.
'normalize.py' is the module that can apply transformations calculated in 'transformation.py' to an image.
'validate.py' is the module for validating radiometric normalization.

Each of the modules referenced above are intended to run on a single band at a time (represented as a numpy array). Examples of more complete functions that can handle image reading and multiple bands are within 'radiometric_normalization/wrappers'. These will be explained below.

'tests' contain unit tests for the functions in the library.

Example Usage

The example below demonstrates the generation of per-band linear transformations that will normalize a candidate image to the mean values of a set of reference images. The candidate and reference images must be 16-bit. Additionally, all of the reference images in the set must have the same number and order bands and pixel dimensions as the candidate image.

candidate_path is a string specifying the location of the candidate image on disk. reference_paths is a list of strings, each specifying the location of a reference image on disk. transformations is a list of tuples, each specifying the gain (first entry) and offset (second entry) that will normalize the respective band of the candidate image.

Below is an example using two Landsat8 tiles: LC08_L1TP_044034_20170427_20170428_01_T1 and LC08_L1TP_044034_20170105_20170218_01_RT.

from radiometric_normalization.wrappers import pif_wrapper
from radiometric_normalization.wrappers import transformation_wrapper
from radiometric_normalization.wrappers import normalize_wrapper
from radiometric_normalization import gimage
from radiometric_normalization import pif

## OPTIONAL
import logging
import numpy
import subprocess
from osgeo import gdal
from radiometric_normalization.wrappers import display_wrapper

logging.basicConfig(level=logging.DEBUG)
##

## OPTIONAL - Cut dataset to colocated sub scenes and create and BGRN image
# LC08_L1TP_044034_20170105_20170218_01_T1 is the older scene and so it is set as the reference.

band_mapping = [{'name': 'blue', 'L8': 'B2'}, {'name': 'green', 'L8':'B3'}, {'name': 'red', 'L8': 'B4'}, {'name': 'nir', 'L8': 'B5'}]

full_candidate_basename = 'LC08_L1TP_044034_20170427_20170428_01_RT'
full_reference_basename = 'LC08_L1TP_044034_20170105_20170218_01_T1'
candidate_basename = 'candidate'
reference_basename = 'reference'
full_candidate_filenames = ['{}_{}.TIF'.format(full_candidate_basename, b['L8']) for b in band_mapping]
candidate_filenames = ['{}_{}.TIF'.format(candidate_basename, b['name']) for b in band_mapping]
full_reference_filenames = ['{}_{}.TIF'.format(full_reference_basename, b['L8']) for b in band_mapping]
reference_filenames = ['{}_{}.TIF'.format(reference_basename, b['name']) for b in band_mapping]

for full_filename, cropped_filename in zip(full_candidate_filenames, candidate_filenames):
    subprocess.check_call(['gdal_translate', '-projwin', '545000', '4136000', '601000', '4084000', full_filename, cropped_filename])

for full_filename, cropped_filename in zip(full_reference_filenames, reference_filenames):
    subprocess.check_call(['gdal_translate', '-projwin', '545000', '4136000', '601000', '4084000', full_filename, cropped_filename])

band_gimgs = {}
for cropped_filename in candidate_filenames:
    band = cropped_filename.split('_')[1].split('.TIF')[0]
    band_gimgs[band] = gimage.load(cropped_filename)

candidate_path = 'candidate.tif'
combined_alpha = numpy.logical_and.reduce([b.alpha for b in band_gimgs.values()])
temporary_gimg = gimage.GImage([band_gimgs[b].bands[0] for b in ['blue', 'green', 'red', 'nir']], combined_alpha, band_gimgs['blue'].metadata)
gimage.save(temporary_gimg, candidate_path)

band_gimgs = {}
for cropped_filename in reference_filenames:
    band = cropped_filename.split('_')[1].split('.TIF')[0]
    band_gimgs[band] = gimage.load(cropped_filename)

reference_path = 'reference.tif'
combined_alpha = numpy.logical_and.reduce([b.alpha for b in band_gimgs.values()])
temporary_gimg = gimage.GImage([band_gimgs[b].bands[0] for b in ['blue', 'green', 'red', 'nir']], combined_alpha, band_gimgs['blue'].metadata)
gimage.save(temporary_gimg, reference_path)
##

parameters = pif.pca_options(threshold=100)
pif_mask = pif_wrapper.generate(candidate_path, reference_path, method='filter_PCA', last_band_alpha=True, method_options=parameters)

## OPTIONAL - Save out the PIF mask
candidate_ds = gdal.Open(candidate_path)
metadata = gimage.read_metadata(candidate_ds)
pif_gimg = gimage.GImage([pif_mask], numpy.ones(pif_mask.shape, dtype=numpy.bool), metadata)
gimage.save(pif_gimg, 'PIF_pixels.tif')
##

transformations = transformation_wrapper.generate(candidate_path, reference_path, pif_mask, method='linear_relationship', last_band_alpha=True)

## OPTIONAL - View the transformations
print transformations
##

normalised_gimg = normalize_wrapper.generate(candidate_path, transformations, last_band_alpha=True)
result_path = 'normalized.tif'
gimage.save(normalised_gimg, result_path)

## OPTIONAL - View the effect on the pixels (SLOW)
from radiometric_normalization.wrappers import display_wrapper
display_wrapper.create_pixel_plots(candidate_path, reference_path, 'Original', limits=[0, 30000], last_band_alpha=True)
display_wrapper.create_pixel_plots(result_path, reference_path, 'Transformed', limits=[0, 30000], last_band_alpha=True)
display_wrapper.create_all_bands_histograms(candidate_path, reference_path, 'Original', x_limits=[4000, 25000], last_band_alpha=True)
display_wrapper.create_all_bands_histograms(result_path, reference_path, 'Transformed', x_limits=[4000, 25000], last_band_alpha=True)
##

In the above example, the 'reference.tif', original candidate scene ('candidate.tif') and radiometrically normalized candidate scene ('normalized.tif') all displayed at the same intensity scale:

Reference	Original	Transformed	PIF (red pixels) over Reference

You can compare the per-band DN to DN plots and the histograms of the reference image with the original candidate image and with the normalized candidate image: