WorldView Mosaic

In [ ]:

# This file demonstrates how to preprocess Worldview3 imagery into a mosaic using spectralmatch.
# Starting from two overlapping Worldview3 images in reflectance, the process includes joint coregistration, global matching, local matching, starting from saved block maps (optional for demonstration purposes), generating seamlines, and marging images, and before vs after statistics.
# This script is set up to perform matching on all .tif files from a folder within the working directory called "Input" e.g. working_directory/Input/*.tif. The easiest way to process your own imagery is to move it inside that folder or change the working_directory to another folder with this structure, alternatively, you can pass in custom lists of image paths.

# This file demonstrates how to preprocess Worldview3 imagery into a mosaic using spectralmatch. # Starting from two overlapping Worldview3 images in reflectance, the process includes joint coregistration, global matching, local matching, starting from saved block maps (optional for demonstration purposes), generating seamlines, and marging images, and before vs after statistics. # This script is set up to perform matching on all .tif files from a folder within the working directory called "Input" e.g. working_directory/Input/*.tif. The easiest way to process your own imagery is to move it inside that folder or change the working_directory to another folder with this structure, alternatively, you can pass in custom lists of image paths.

In [ ]:

import os
from spectralmatch import (
    align_rasters,
    compare_before_after_all_images,
    compare_image_spectral_profiles_pairs,
    compare_spatial_spectral_difference_band_average,
    joint_coregistration,
    global_regression,
    local_block_adjustment,
    mask_rasters,
    merge_rasters,
    search_paths,
    voronoi_center_seamline,
    weighted_seamline,
)

# Important: If this does not automatically find the correct CWD, manually copy the path to the /data_worldview folder
working_directory = os.path.join(os.getcwd(), "data_worldview")
print(working_directory)

input_folder = os.path.join(working_directory, "Input")
coregistration_folder = os.path.join(working_directory, "Coregistration")
global_folder = os.path.join(working_directory, "GlobalMatch")
local_folder = os.path.join(working_directory, "LocalMatch")
clipped_folder = os.path.join(working_directory, "Clipped")
stats_folder = os.path.join(working_directory, "Stats")

window_size = 1024
image_threads = 3 # Dask: None | int
io_threads = 3
tile_threads = 3
debug_mode = True

concurrent_processing_backend = "process_pool"  # process_pool | dask
dask_scheduler = None  # None | ("file", f"/tmp/spectralmatch-dask-1.json") | ("address", "tcp://scheduler:8786")

import os from spectralmatch import ( align_rasters, compare_before_after_all_images, compare_image_spectral_profiles_pairs, compare_spatial_spectral_difference_band_average, joint_coregistration, global_regression, local_block_adjustment, mask_rasters, merge_rasters, search_paths, voronoi_center_seamline, weighted_seamline, ) # Important: If this does not automatically find the correct CWD, manually copy the path to the /data_worldview folder working_directory = os.path.join(os.getcwd(), "data_worldview") print(working_directory) input_folder = os.path.join(working_directory, "Input") coregistration_folder = os.path.join(working_directory, "Coregistration") global_folder = os.path.join(working_directory, "GlobalMatch") local_folder = os.path.join(working_directory, "LocalMatch") clipped_folder = os.path.join(working_directory, "Clipped") stats_folder = os.path.join(working_directory, "Stats") window_size = 1024 image_threads = 3 # Dask: None | int io_threads = 3 tile_threads = 3 debug_mode = True concurrent_processing_backend = "process_pool" # process_pool | dask dask_scheduler = None # None | ("file", f"/tmp/spectralmatch-dask-1.json") | ("address", "tcp://scheduler:8786")

In [ ]:

joint_coregistration(
    input_images=input_folder,
    output_images=coregistration_folder,
    tap=True,
    resolution="highest",
    debug_logs=debug_mode,
    window_size=window_size,
    image_threads=image_threads,
    io_threads=io_threads,
    tile_threads=tile_threads,
    concurrent_processing_backend=concurrent_processing_backend,
    dask_scheduler=dask_scheduler,
)

joint_coregistration( input_images=input_folder, output_images=coregistration_folder, tap=True, resolution="highest", debug_logs=debug_mode, window_size=window_size, image_threads=image_threads, io_threads=io_threads, tile_threads=tile_threads, concurrent_processing_backend=concurrent_processing_backend, dask_scheduler=dask_scheduler, )

In [ ]:

# align_rasters(
#     input_images=input_folder,
#     output_images=coregistration_folder,
#     tap=True,
#     resolution="highest",
#     debug_logs=debug_mode,
#     window_size=window_size,
#     image_threads=image_threads,
#     io_threads=io_threads,
#     tile_threads=tile_threads,
#     concurrent_processing_backend=concurrent_processing_backend,
#     dask_scheduler=dask_scheduler,
# )

# align_rasters( # input_images=input_folder, # output_images=coregistration_folder, # tap=True, # resolution="highest", # debug_logs=debug_mode, # window_size=window_size, # image_threads=image_threads, # io_threads=io_threads, # tile_threads=tile_threads, # concurrent_processing_backend=concurrent_processing_backend, # dask_scheduler=dask_scheduler, # )

In [ ]:

inz_output_path = os.path.join(global_folder, "INZ", "$_to_$_INZ.tif")

global_regression(
    input_images=coregistration_folder, # Automatically searches for all *.tif files if passed this way
    output_images=global_folder,
    debug_logs=debug_mode,
    window_size=window_size,
    image_threads=image_threads,
    io_threads=io_threads,
    tile_threads=tile_threads,
    concurrent_processing_backend=concurrent_processing_backend,
    dask_scheduler=dask_scheduler,
    save_as_cog=True,
    # custom_nodata_value=-9999,
    pif_method="flood_from_match_points",
    estimate_stats=True,
    # specify_model_images=("include", ['Worldview_2016-09-22']), # Global matching all input images to the spectral profile of any number of specified images (regression will still be based on overlapping areas, however, only the *included* images statistics will influence the solution)
    # custom_mean_factor=3, # Default is 1; 3 often works better to 'move' the spectral mean of images closer together (applied when creating model)
    # custom_std_factor=3,
    save_adjustments=os.path.join(
        global_folder, "GlobalAdjustments.json"
    ),  # Start from precomputed statistics for images whole and overlap stats
    pif_save_inz=inz_output_path,  # Saves one INZ raster per overlap pair; first $ = main/sensed image name, second $ = reference image name
    # load_adjustments=os.path.join(global_folder, "GlobalAdjustments.json"), # Load Statistics
)

inz_output_path = os.path.join(global_folder, "INZ", "$_to_$_INZ.tif") global_regression( input_images=coregistration_folder, # Automatically searches for all *.tif files if passed this way output_images=global_folder, debug_logs=debug_mode, window_size=window_size, image_threads=image_threads, io_threads=io_threads, tile_threads=tile_threads, concurrent_processing_backend=concurrent_processing_backend, dask_scheduler=dask_scheduler, save_as_cog=True, # custom_nodata_value=-9999, pif_method="flood_from_match_points", estimate_stats=True, # specify_model_images=("include", ['Worldview_2016-09-22']), # Global matching all input images to the spectral profile of any number of specified images (regression will still be based on overlapping areas, however, only the *included* images statistics will influence the solution) # custom_mean_factor=3, # Default is 1; 3 often works better to 'move' the spectral mean of images closer together (applied when creating model) # custom_std_factor=3, save_adjustments=os.path.join( global_folder, "GlobalAdjustments.json" ), # Start from precomputed statistics for images whole and overlap stats pif_save_inz=inz_output_path, # Saves one INZ raster per overlap pair; first $ = main/sensed image name, second $ = reference image name # load_adjustments=os.path.join(global_folder, "GlobalAdjustments.json"), # Load Statistics )

In [ ]:

reference_map_path = os.path.join(local_folder, "ReferenceBlockMap", "ReferenceBlockMap.tif")
local_maps_path = os.path.join(local_folder, "LocalBlockMap", "$_LocalBlockMap.tif")
# Only search for saved maps when enabling load_block_maps below.
# searched_paths = search_paths(os.path.join(local_folder, "LocalBlockMap", "*.tif"))

local_block_adjustment(
    input_images=global_folder,
    output_images=local_folder,
    debug_logs=debug_mode,
    window_size=window_size,
    image_threads=image_threads,
    io_threads=io_threads,
    tile_threads=tile_threads,
    concurrent_processing_backend=concurrent_processing_backend,
    dask_scheduler=dask_scheduler,
    save_as_cog=True,
    # custom_nodata_value=-9999,
    correction_method="offset",
    number_of_blocks=50,  # Target number of blocks
    # override_bounds_canvas_coords = (193011.1444011169369332, 2184419.3597142999060452, 205679.2836037494416814, 2198309.8632259583100677), # Local match with a larger canvas than images bounds (perhaps to anticipate adding additional imagery so you don't have to recalculate local block maps each rematch)
    save_block_maps=(reference_map_path, local_maps_path),
    # load_block_maps=(reference_map_path, searched_paths), # Local match from saved block maps (this code just passes in local maps, but if a reference map is passed in, it will match images to the reference map without recomputing it)
)

reference_map_path = os.path.join(local_folder, "ReferenceBlockMap", "ReferenceBlockMap.tif") local_maps_path = os.path.join(local_folder, "LocalBlockMap", "$_LocalBlockMap.tif") # Only search for saved maps when enabling load_block_maps below. # searched_paths = search_paths(os.path.join(local_folder, "LocalBlockMap", "*.tif")) local_block_adjustment( input_images=global_folder, output_images=local_folder, debug_logs=debug_mode, window_size=window_size, image_threads=image_threads, io_threads=io_threads, tile_threads=tile_threads, concurrent_processing_backend=concurrent_processing_backend, dask_scheduler=dask_scheduler, save_as_cog=True, # custom_nodata_value=-9999, correction_method="offset", number_of_blocks=50, # Target number of blocks # override_bounds_canvas_coords = (193011.1444011169369332, 2184419.3597142999060452, 205679.2836037494416814, 2198309.8632259583100677), # Local match with a larger canvas than images bounds (perhaps to anticipate adding additional imagery so you don't have to recalculate local block maps each rematch) save_block_maps=(reference_map_path, local_maps_path), # load_block_maps=(reference_map_path, searched_paths), # Local match from saved block maps (this code just passes in local maps, but if a reference map is passed in, it will match images to the reference map without recomputing it) )

In [ ]:

# Option 1: Voronoi center seamlines from the raster footprints

voronoi_center_seamline(
    input_images=local_folder,
    output_mask=os.path.join(working_directory, "ImageMasks.gpkg"),
    image_field_name="image",
    debug_logs=debug_mode,
    debug_vectors_path=os.path.join(working_directory, "DebugVectors.gpkg"),
)

# Option 2: Weighted seamlines from a polygon layer that already contains ranking attributes. This is useful when another workflow prepares footprint polygons plus metadata fields for ranking. Rank placeholders support direct fields like {quality_score}.
# weighted_seamline(
#     input_polygons=os.path.join(working_directory, "SeamlineMetadata.gpkg"),
#     output_mask=os.path.join(working_directory, "ImageMasks.gpkg"),
#     input_layer="footprints",
#     image_field_name="image",
#     rank_function="{quality_score} - {cloud_cover}",
#     rank_descending=True,
#     debug_logs=debug_mode,
# )

# Option 1: Voronoi center seamlines from the raster footprints voronoi_center_seamline( input_images=local_folder, output_mask=os.path.join(working_directory, "ImageMasks.gpkg"), image_field_name="image", debug_logs=debug_mode, debug_vectors_path=os.path.join(working_directory, "DebugVectors.gpkg"), ) # Option 2: Weighted seamlines from a polygon layer that already contains ranking attributes. This is useful when another workflow prepares footprint polygons plus metadata fields for ranking. Rank placeholders support direct fields like {quality_score}. # weighted_seamline( # input_polygons=os.path.join(working_directory, "SeamlineMetadata.gpkg"), # output_mask=os.path.join(working_directory, "ImageMasks.gpkg"), # input_layer="footprints", # image_field_name="image", # rank_function="{quality_score} - {cloud_cover}", # rank_descending=True, # debug_logs=debug_mode, # )

In [ ]:

mask_rasters(
    input_images=local_folder,
    output_images=clipped_folder,
    vector_mask=("include", os.path.join(working_directory, "ImageMasks.gpkg"), "image"),
    debug_logs=debug_mode,
    window_size=window_size,
    image_threads=image_threads,
    io_threads=io_threads,
    tile_threads=tile_threads,
    concurrent_processing_backend=concurrent_processing_backend,
    dask_scheduler=dask_scheduler,
)

mask_rasters( input_images=local_folder, output_images=clipped_folder, vector_mask=("include", os.path.join(working_directory, "ImageMasks.gpkg"), "image"), debug_logs=debug_mode, window_size=window_size, image_threads=image_threads, io_threads=io_threads, tile_threads=tile_threads, concurrent_processing_backend=concurrent_processing_backend, dask_scheduler=dask_scheduler, )

In [ ]:

merge_rasters(
    input_images=clipped_folder,
    # output_image_path=os.path.join(working_directory, "MergedImage.tif"), # Use this for single tif output and set output_tiles=False
    output_image_path=os.path.join(working_directory, "MergedImage"),
    debug_logs=debug_mode,
    window_size=window_size,
    io_threads=io_threads,
    tile_threads=tile_threads,
    build_overviews=True,
    output_tiles=True
)

merge_rasters( input_images=clipped_folder, # output_image_path=os.path.join(working_directory, "MergedImage.tif"), # Use this for single tif output and set output_tiles=False output_image_path=os.path.join(working_directory, "MergedImage"), debug_logs=debug_mode, window_size=window_size, io_threads=io_threads, tile_threads=tile_threads, build_overviews=True, output_tiles=True )

In [ ]:

# Compare image spectral profiles pairs
image_pairs = {
    os.path.splitext(os.path.basename(b))[0]: [b, a]
    for b, a in zip(search_paths(os.path.join(input_folder, "*.tif")), search_paths(os.path.join(local_folder, "*.tif")))
}

compare_image_spectral_profiles_pairs(
    image_pairs,
    os.path.join(stats_folder, "ImageSpectralProfilesPairs.png"),
    title="Comparison of Image Spectral Profile Pairs",
    xlabel="Band",
    ylabel="Reflectance",
)

# Compare spatial spectral difference band average
before_paths, after_paths = zip(*zip(search_paths(os.path.join(input_folder, "*.tif")), search_paths(os.path.join(local_folder, "*.tif"))))

for before_path, after_path in zip(before_paths, after_paths):
    compare_spatial_spectral_difference_band_average(
        input_images=[before_path, after_path],
        output_figure_path=os.path.join(
            stats_folder,
            f"PixelChange_{os.path.splitext(os.path.basename(before_path))[0]}.png"
        ),
        title="Input to Output Comparison of Pixel Change",
        diff_label="Pixel Difference",
        subtitle=f"Image: {os.path.splitext(os.path.basename(before_path))[0]}",
    )

# Compare before after all images
compare_before_after_all_images(
    input_images_1=search_paths(os.path.join(input_folder, "*.tif")),
    input_images_2=search_paths(os.path.join(local_folder, "*.tif")),
    output_figure_path=os.path.join(stats_folder, "CompareBeforeAfterAllImages.png"),
    image_names=[os.path.splitext(os.path.basename(p))[0] for p in search_paths(os.path.join(input_folder, "*.tif"))],
    title="Comparison of Before to After of all Images",
    ylabel_1="Before",
    ylabel_2="After",
)

# Compare image spectral profiles pairs image_pairs = { os.path.splitext(os.path.basename(b))[0]: [b, a] for b, a in zip(search_paths(os.path.join(input_folder, "*.tif")), search_paths(os.path.join(local_folder, "*.tif"))) } compare_image_spectral_profiles_pairs( image_pairs, os.path.join(stats_folder, "ImageSpectralProfilesPairs.png"), title="Comparison of Image Spectral Profile Pairs", xlabel="Band", ylabel="Reflectance", ) # Compare spatial spectral difference band average before_paths, after_paths = zip(*zip(search_paths(os.path.join(input_folder, "*.tif")), search_paths(os.path.join(local_folder, "*.tif")))) for before_path, after_path in zip(before_paths, after_paths): compare_spatial_spectral_difference_band_average( input_images=[before_path, after_path], output_figure_path=os.path.join( stats_folder, f"PixelChange_{os.path.splitext(os.path.basename(before_path))[0]}.png" ), title="Input to Output Comparison of Pixel Change", diff_label="Pixel Difference", subtitle=f"Image: {os.path.splitext(os.path.basename(before_path))[0]}", ) # Compare before after all images compare_before_after_all_images( input_images_1=search_paths(os.path.join(input_folder, "*.tif")), input_images_2=search_paths(os.path.join(local_folder, "*.tif")), output_figure_path=os.path.join(stats_folder, "CompareBeforeAfterAllImages.png"), image_names=[os.path.splitext(os.path.basename(p))[0] for p in search_paths(os.path.join(input_folder, "*.tif"))], title="Comparison of Before to After of all Images", ylabel_1="Before", ylabel_2="After", )