Note

This page was generated from a Jupyter notebook. To run and interact with it, you can download it here.

Basic correction and evaluation

This example corrects a small set of aerial images and evaluates the results.

Images are taken from the homonim test data: Aerials were obtained from NGI, and subsequently downsampled and re-scaled to uint8 precision to reduce size. Satellite images were acquired with geedim, and re-scaled to uint8 precision.

The notebook uses the APIs of homonim and other packages. CLI commands equivalent to the API code snippets are given in the comments where possible.

Setup

gdal and matplotlib are required to run the notebook. You can uncomment the cell below to install them, if they aren’t installed already.

# matplotlib and gdal should be installed if they aren't already.
# import sys
# if 'conda' in sys.prefix:
#     # install into the conda environment the notebook is being run from
#     !conda install --yes --prefix {sys.prefix} -c conda-forge gdal matplotlib
# else:
#     # install into the python environment the notebook is being run from
#     !{sys.executable} -m pip install gdal matplotlib

import logging
from pathlib import Path
from matplotlib import pyplot
import rasterio as rio
from rasterio.plot import show
import numpy as np
from homonim import RasterFuse, RasterCompare, Model

logging.basicConfig(level=logging.ERROR)

# create URLS of homonim source and reference images.
src_paths = [
    f'https://raw.githubusercontent.com/leftfield-geospatial/homonim/main/'
    f'tests/data/source/ngi_rgb_byte_{i}.tif'
    for i in range(1, 5)
]
ref_path = (
    'https://raw.githubusercontent.com/leftfield-geospatial/homonim/main/'
    'tests/data/reference/sentinel2_b432_byte.tif'
)

Surface reflectance correction

Here we fuse the aerial images with a Sentinel-2 reference, using the RasterFuse class. Source and reference images are RGB, so bands are matched automatically without specifying the src_bands or ref_bands parameters. Fusion uses the gain-blk-offset model, with a kernel shape of 5x5 pixels. These standard settings work reasonably well in a variety of contexts.

# corrected file names corresponding to source names
corr_paths = [f'corr_{Path(src_path).stem}.tif' for src_path in src_paths]

for src_path, corr_path in zip(src_paths, corr_paths):
    with RasterFuse(src_path, ref_path) as raster_fuse:
        raster_fuse.process(
            corr_path, Model.gain_blk_offset, (5, 5), overwrite=True
        )

    # equivalent homonim command line:
    # !homonim fuse -m gain-blk-offset -k 5 5 -o {src_path} {ref_path}

Visualisation

Next, we create a VRT mosaic of the corrected images that will be used here, and in the evaluation step.

# strictly, one should avoid using GDAL and rasterio together, but it doesn't
# create conflicts here
from osgeo import gdal

corr_mosaic_path = 'corr_ngi_mosaic_rgb_byte.vrt'
ds = gdal.BuildVRT(corr_mosaic_path, corr_paths)
ds.FlushCache()

Now we can display a zoomed-in area of the source mosaic, reference image and corrected mosaic.

# VRT mosaic of the source files
src_mosaic_path = (
    'https://raw.githubusercontent.com/leftfield-geospatial/homonim/main/'
    'tests/data/source/ngi_mosaic_rgb_byte.vrt'
)

fig, axes = pyplot.subplots(
    1, 3, sharex=True, sharey=True, tight_layout=True, figsize=(9, 4), dpi=92
)

# loop over the source, reference and corrected images etc
for im_file, scale, ax, label in zip(
    [src_mosaic_path, ref_path, corr_mosaic_path],
    [255, 150, 150],
    axes,
    ['Source', 'Reference', 'Corrected'],
):
    # read, scale and display the image
    with rio.open(im_file, 'r') as ds:
        array = ds.read(out_dtype='float32') / scale

        show(array, transform=ds.transform, ax=ax, interpolation='nearest')

    ax.set_title(label)
    ax.set_xlim(-5.75e4, -5.50e4)  # zoom in
    ax.set_ylim(-3.733e6, -3.730e6)
    ax.axis('off')

# fig.savefig('../case_studies/basic_correction-src_ref_corr.jpg', dpi=92)

The corrected mosaic is free of seamlines and similar in colour to the reference.

Evaluation

In this last section, we compare the source and corrected similarity with another reference. The comparison gives an indication of the change in surface reflectance accuracy due to correction. Rather than compare with the Sentinel-2 reference used for correction, we compare with an independent Landsat-8 reference.

RasterCompare is used to produce tables of comparison statistics.

# URL of the Landsat-8 reference
cmp_ref_path = (
    'https://raw.githubusercontent.com/leftfield-geospatial/homonim/main/'
    'tests/data/reference/landsat8_byte.tif'
)
print(RasterCompare.schema_table())

# loop over the source and corrected image files
for im_path, im_label in zip(
    [src_mosaic_path, corr_mosaic_path],
    ['Source', 'Corrected'],
):
    with RasterCompare(im_path, cmp_ref_path) as compare:
        # print a table of comparison statistics (the typical way of using
        # RasterCompare)
        stats_dict = compare.process()
        print(f'{im_label} comparison:\n\n' + compare.stats_table(stats_dict))

    # equivalent homonim command line:
    # !homonim compare {im_path} {cmp_ref_path}

ABBREV   DESCRIPTION
-------- -----------------------------------------
r²       Pearson's correlation coefficient squared
RMSE     Root Mean Square Error
rRMSE    Relative RMSE (RMSE/mean(ref))
N        Number of pixels

Source comparison:

  Band    r²   RMSE   rRMSE     N
------ ----- ------ ------- -----
 SR_B4 0.639 79.247   1.746 76143
 SR_B3 0.549 83.168   1.859 76143
 SR_B2 0.409 94.104   3.169 76143
  Mean 0.533 85.506   2.258 76143

Corrected comparison:

  Band    r²   RMSE   rRMSE     N
------ ----- ------ ------- -----
 SR_B4 0.938  9.815   0.216 76143
 SR_B3 0.922  6.892   0.154 76143
 SR_B2 0.862 30.966   1.043 76143
  Mean 0.907 15.891   0.471 76143

The tables (and \(r^2\) values in particular) show that correction has improved similarity with the Landsat-8 reference.