Statistical Figures

compare_before_after_all_images(input_images_1, input_images_2, output_figure_path, title, ylabel_1, ylabel_2, image_names=None)

Creates a two-row image grid to compare before-and-after raster pairs with consistent per-row contrast stretching. Each column shows a pair of aligned images with transparent nodata. Supports 1- and 3-band rasters.

Parameters:

Name	Type	Description	Default
input_images_1	list	List of file paths to the "before" images (top row).	required
input_images_2	list	List of file paths to the "after" images (bottom row).	required
output_figure_path	str	Destination path to save the output PNG figure.	required
title	str	Title of the entire figure.	required
ylabel_1	str	Y-axis label for the top row.	required
ylabel_2	str	Y-axis label for the bottom row.	required
image_names	list	List of image names to use as column titles. Must match the number of image pairs.	None

Raises:

Type	Description
AssertionError	If input lists have mismatched lengths or if image_names does not match image count.

Output

Saves a PNG file with the comparison figure.

Source code in spectralmatch/statistics.py

def compare_before_after_all_images(
    input_images_1: list,
    input_images_2: list,
    output_figure_path: str,
    title: str,
    ylabel_1: str,
    ylabel_2: str,
    image_names: list = None,
):
    """
    Creates a two-row image grid to compare before-and-after raster pairs with consistent per-row contrast stretching. Each column shows a pair of aligned images with transparent nodata. Supports 1- and 3-band rasters.

    Args:
        input_images_1 (list): List of file paths to the "before" images (top row).
        input_images_2 (list): List of file paths to the "after" images (bottom row).
        output_figure_path (str): Destination path to save the output PNG figure.
        title (str): Title of the entire figure.
        ylabel_1 (str): Y-axis label for the top row.
        ylabel_2 (str): Y-axis label for the bottom row.
        image_names (list, optional): List of image names to use as column titles. Must match the number of image pairs.

    Raises:
        AssertionError: If input lists have mismatched lengths or if `image_names` does not match image count.

    Output:
        Saves a PNG file with the comparison figure.
    """
    _print_step_start("compare_before_after_all_images")
    def read_as_3band(ds):
        count = ds.RasterCount
        if count >= 3:
            arr = np.stack(
                [ds.GetRasterBand(i + 1).ReadAsArray() for i in range(3)]
            )
        else:
            single = ds.GetRasterBand(1).ReadAsArray()
            arr = np.repeat(single[np.newaxis, ...], 3, axis=0)
        return arr.astype("float32")

    def compute_row_stretch(paths):
        all_valid = [[] for _ in range(3)]
        for path in paths:
            _print_image_start(path, output_figure_path)
            ds = gdal.Open(path, gdal.GA_ReadOnly)
            if ds is None:
                raise RuntimeError(f"Failed to open {path}")

            nodata = ds.GetRasterBand(1).GetNoDataValue()
            img = read_as_3band(ds)

            mask = np.full(img.shape[1:], False)
            if nodata is not None:
                for b in range(img.shape[0]):
                    mask |= img[b] == nodata

            for b in range(img.shape[0]):
                all_valid[b].append(img[b][~mask])

            ds = None

        return [
            np.percentile(np.concatenate(valid), (2, 98))
            if valid else (0, 1)
            for valid in all_valid
        ]

    assert len(input_images_1) == len(input_images_2)
    if image_names:
        assert len(image_names) == len(input_images_1)

    os.makedirs(os.path.dirname(output_figure_path), exist_ok=True)

    num_images = len(input_images_1)
    fig = plt.figure(figsize=(5 * num_images, 10))
    gs = gridspec.GridSpec(2, num_images + 1, width_ratios=[0.05] + [1] * num_images)

    stretch_1 = compute_row_stretch(input_images_1)
    stretch_2 = compute_row_stretch(input_images_2)

    for col_idx, (path1, path2) in enumerate(zip(input_images_1, input_images_2)):
        for row_idx, (path, stretch) in enumerate(
            [(path1, stretch_1), (path2, stretch_2)]
        ):
            ax = fig.add_subplot(gs[row_idx, col_idx + 1])

            ds = gdal.Open(path, gdal.GA_ReadOnly)
            if ds is None:
                raise RuntimeError(f"Failed to open {path}")

            nodata = ds.GetRasterBand(1).GetNoDataValue()
            img = read_as_3band(ds)

            mask = np.full(img.shape[1:], False)
            if nodata is not None:
                for b in range(img.shape[0]):
                    mask |= img[b] == nodata

            for b in range(img.shape[0]):
                vmin, vmax = stretch[b]
                img[b] = np.clip((img[b] - vmin) / (vmax - vmin), 0, 1)

            img = img.transpose(1, 2, 0)
            alpha = (~mask).astype("float32")
            rgba = np.dstack((img, alpha))

            ax.imshow(rgba)

            if row_idx == 0 and image_names:
                ax.set_title(image_names[col_idx])

            ax.axis("off")
            ds = None

    for i, label in enumerate([ylabel_1, ylabel_2]):
        ax = fig.add_subplot(gs[i, 0])
        ax.set_ylabel(label, fontsize=12, rotation=90, labelpad=10, va="center")
        ax.tick_params(left=False, labelleft=False)
        ax.set_xticks([])
        ax.set_yticks([])
        ax.set_frame_on(False)

    fig.suptitle(title, fontsize=16)
    plt.tight_layout(rect=[0, 0, 1, 0.95])
    plt.savefig(output_figure_path, dpi=300)
    plt.close()

    print(f"Saved: {os.path.splitext(os.path.basename(output_figure_path))[0]}")
    for completed, path in enumerate([*input_images_1, *input_images_2], 1):
        _print_image_completed(path, completed, len(input_images_1) + len(input_images_2))

compare_image_spectral_profiles_pairs(image_groups_dict, output_figure_path, title, xlabel, ylabel, line_width=1, estimate_stats=True)

Plots paired spectral profiles for before-and-after image comparisons.

Parameters:

Name	Type	Description	Default
image_groups_dict	dict	Mapping of labels to image path pairs (before, after): {'Image A': [ '/image/before/a.tif', 'image/after/a.tif' ], 'Image B': [ '/image/before/b.tif', '/image/after/b.tif' ]}	required
output_figure_path	str	Path to save the resulting comparison figure.	required
title	str	Title of the plot.	required
xlabel	str	X-axis label.	required
ylabel	str	Y-axis label.	required
line_width	float	Width of the spectral profiles lines. Default is 1.	1
estimate_stats	bool	Whether to estimate band statistics. Default is True.	True

Outputs

Saves a spectral comparison plot showing pre- and post-processing profiles.

Source code in spectralmatch/statistics.py

def compare_image_spectral_profiles_pairs(
    image_groups_dict: dict,
    output_figure_path: str,
    title: str,
    xlabel: str,
    ylabel: str,
    line_width: float = 1,
    estimate_stats: bool = True,
):
    """
    Plots paired spectral profiles for before-and-after image comparisons.

    Args:
        image_groups_dict (dict): Mapping of labels to image path pairs (before, after):
            {'Image A': [
                '/image/before/a.tif',
                'image/after/a.tif'
            ],
            'Image B': [
                '/image/before/b.tif',
                '/image/after/b.tif'
            ]}
        output_figure_path (str): Path to save the resulting comparison figure.
        title (str): Title of the plot.
        xlabel (str): X-axis label.
        ylabel (str): Y-axis label.
        line_width (float, optional): Width of the spectral profiles lines. Default is 1.
        estimate_stats (bool, optional): Whether to estimate band statistics. Default is True.

    Outputs:
        Saves a spectral comparison plot showing pre- and post-processing profiles.
    """

    _print_step_start("compare_image_spectral_profiles_pairs")
    os.makedirs(os.path.dirname(output_figure_path), exist_ok=True)
    plt.figure(figsize=(10, 6))
    colors = itertools.cycle(plt.cm.tab10.colors)
    processed_paths = []

    for label, group in image_groups_dict.items():
        if len(group) == 2:
            image_path1, image_path2 = group
            color = next(colors)

            for i, image_path in enumerate([image_path1, image_path2]):
                _print_image_start(image_path, output_figure_path)
                ds = gdal.Open(image_path, gdal.GA_ReadOnly)
                if ds is None:
                    continue
                processed_paths.append(image_path)

                num_bands = ds.RasterCount
                mean_spectral = np.zeros(num_bands, dtype=float)

                for b in range(1, num_bands + 1):
                    band = ds.GetRasterBand(b)

                    # (min, max, mean, std)
                    stats = band.GetStatistics(False, estimate_stats)

                    # stats[2] is the mean
                    mean_spectral[b - 1] = stats[2]

                ds = None

                bands = np.arange(1, num_bands + 1)
                linestyle = "dashed" if i == 0 else "solid"

                plt.plot(
                    bands,
                    mean_spectral,
                    linestyle=linestyle,
                    color=color,
                    linewidth=line_width,
                    label=f"{label} - {'Before' if i == 0 else 'After'}",
                )

                ds = None

    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)
    plt.legend(
        loc="upper center",
        bbox_to_anchor=(0.5, -0.15),
        ncol=2,
        frameon=True,
    )
    plt.tight_layout()
    plt.grid(True)
    plt.xticks(np.arange(1, num_bands + 1, 1))
    plt.savefig(output_figure_path, dpi=300)
    plt.close()
    print(f"Saved: {os.path.splitext(os.path.basename(output_figure_path))[0]}")
    for completed, path in enumerate(processed_paths, 1):
        _print_image_completed(path, completed, len(processed_paths))

compare_spatial_spectral_difference_band_average(input_images, output_figure_path, title, diff_label, subtitle, scale=None, *, resampling_method='bilinear')

Reprojects the after image onto the before image's geographic grid and plots the signed mean spectral difference (after minus before) across matching bands.

Parameters:

Name	Type	Description	Default
input_images	list	Two georeferenced raster paths [before, after] with matching band counts; dimensions, extents, pixel sizes, and CRSs may differ. The before image defines the comparison CRS, extent, and pixel grid; pixels must be valid in both images across all bands to contribute.	required
output_figure_path	str	Path to save the resulting difference image (PNG).	required
title	str	Title for the plot.	required
diff_label	str	Label for the colorbar.	required
subtitle	str	Subtitle text shown below the image.	required
scale	tuple	Tuple (vmin, vmax) to fix the color scale. Centered at 0.	None
resampling_method	Literal['nearest', 'bilinear', 'cubic', 'lanczos']	Resampling used to project the after image onto the before grid: nearest, bilinear (default), cubic, or lanczos. Areas outside shared valid coverage remain transparent.	'bilinear'

Raises:

Type	Description
ValueError	If the input count or resampling method is invalid, band counts differ, georeferencing is missing, or the images have no shared valid pixels.

Source code in spectralmatch/statistics.py

def compare_spatial_spectral_difference_band_average(
    input_images: list,
    output_figure_path: str,
    title: str,
    diff_label: str,
    subtitle: str,
    scale: tuple = None,
    *,
    resampling_method: Literal["nearest", "bilinear", "cubic", "lanczos"] = "bilinear",
):
    """
    Reprojects the after image onto the before image's geographic grid and plots the signed mean spectral difference (after minus before) across matching bands.

    Args:
        input_images (list): Two georeferenced raster paths [before, after] with matching band counts; dimensions, extents, pixel sizes, and CRSs may differ. The before image defines the comparison CRS, extent, and pixel grid; pixels must be valid in both images across all bands to contribute.
        output_figure_path (str): Path to save the resulting difference image (PNG).
        title (str): Title for the plot.
        diff_label (str): Label for the colorbar.
        subtitle (str): Subtitle text shown below the image.
        scale (tuple, optional): Tuple (vmin, vmax) to fix the color scale. Centered at 0.
        resampling_method: Resampling used to project the after image onto the before grid: nearest, bilinear (default), cubic, or lanczos. Areas outside shared valid coverage remain transparent.

    Raises:
        ValueError: If the input count or resampling method is invalid, band counts differ, georeferencing is missing, or the images have no shared valid pixels.
    """
    _print_step_start("compare_spatial_spectral_difference_band_average")
    if len(input_images) != 2:
        raise ValueError("input_images must be a list of exactly two image paths.")
    if resampling_method not in {"nearest", "bilinear", "cubic", "lanczos"}:
        raise ValueError("resampling_method must be nearest, bilinear, cubic, or lanczos.")

    path1, path2 = input_images
    _print_image_start(input_images, output_figure_path)

    ds1 = gdal.Open(path1, gdal.GA_ReadOnly)
    ds2 = gdal.Open(path2, gdal.GA_ReadOnly)

    if ds1 is None or ds2 is None:
        raise RuntimeError("Failed to open one or both input images.")

    try:
        mean_diff = _projected_mean_spectral_difference(ds1, ds2, resampling_method)
    finally:
        ds1 = None
        ds2 = None

    os.makedirs(os.path.dirname(os.path.abspath(output_figure_path)), exist_ok=True)
    fig, ax = plt.subplots(figsize=(10, 6), constrained_layout=True)

    vmin, vmax = scale if scale else (np.nanmin(mean_diff), np.nanmax(mean_diff))
    max_abs = max(abs(vmin), abs(vmax))
    im = ax.imshow(mean_diff, cmap="coolwarm", vmin=-max_abs, vmax=max_abs)

    cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    cbar.set_label(diff_label)

    ax.set_title(title, fontsize=14, pad=12)
    if subtitle:
        ax.text(0.5, -0.1, subtitle, fontsize=10, ha="center", transform=ax.transAxes)

    ax.axis("off")
    plt.savefig(output_figure_path, dpi=300, bbox_inches="tight")
    plt.close()

    print(f"Saved: {os.path.splitext(os.path.basename(output_figure_path))[0]}")
    _print_image_completed(input_images, 1, 1)