Missing Areas After Automatic Georeferencing and Bounding Box Issue

[LiamDiane]

• arosics version:
• Python version:
• Operating System:

After automatic georeferencing, parts of the corrected image that extend beyond the original image's bounding box are missing. The bounding box of the corrected image is not adjusted as expected and instead matches the bounding box of the target image. As a result, areas that extend beyond the target image's bounding box are not visible. How can this be addressed?

[danschef]

Could you please post the version and the exact call of AROSICS you are using?

[LiamDiane]

The version I am using was downloaded from arosics-1.12.1-py311h1ea47a8_0.conda. Below is the output from AROSICS. Thank you.

Automatically detected nodata value for GeoArray_CoReg 'WD_HSI_卫图_Level_15': 0.0
Calculating footprint polygon and actual data corner coordinates for reference image...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:08
C:\Users\lyc\anaconda3\envs\qgis\Lib\site-packages\arosics\CoReg.py:123: UserWarning: The footprint of the reference image contains multiple separate image parts. AROSICS will only process the largest image part.
  warnings.warn(f'The footprint of the {self.imName} contains multiple separate image parts. '
Bounding box of calculated footprint for reference image:
        (3232715.84762773, 4903849.910159007, 3330892.278822106, 5003784.799127045)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Bounding box of calculated footprint for image to be shifted:
        (3237740.123953076, 4904274.882939435, 3333595.4786503827, 5001916.962132358)
Matching window position (X,Y): 3284116.6994189164/4953494.37770168
Target window size (2048, 2048) not possible due to too small overlap area or window position too close to an image edge. New matching window size: (1290, 1290).
Detected integer shifts (X/Y):                            -80/0
Detected subpixel shifts (X/Y):                           -0.33208298683166504/0.35782933235168457
Calculated total shifts in fft pixel units (X/Y):         -80.33208298683167/0.35782933235168457
Calculated total shifts in reference pixel units (X/Y):   -80.33208298683167/0.35782933235168457
Calculated total shifts in target pixel units (X/Y):      -80.33208298683167/0.35782933235168457
Calculated map shifts (X,Y):                              -3879.224336663261/-17.27952522598207
Calculated absolute shift vector length in map units:     3879.2628212268055
Calculated angle of shift vector in degrees from North:   89.74478472147447
Original map info: ['Mercator_1SP', 1.0, 1.0, 3237595.25439938, 5001965.25198359, 48.2898512329, 48.2898512329]
Updated map info:  ['Mercator_1SP', 1.0, 1.0, '3233716.030062717', '5001947.972458364', 48.2898512329, 48.2898512329]
Image similarity within the matching window (SSIM before/after correction): 0.4880 => 0.3235
Estimated reliability of the calculated shifts:  79.5%
Correcting geometric shifts...
Warping progress     |==================================================| 100.0% Complete  => 0:00:11
Writing GeoArray of size (2023, 1987, 90) to C:/Users/lyc/OneDrive/Desktop/georeference-output/test_auto.tif.
success

[LiamDiane]

I resolved the issue temporarily by expanding the bounding box of the image to be corrected using warp before performing automatic correction with arosics. Below is my code:

# Input image paths
im_ref = "D:/LASAC/georeference/WD_HSI_卫图/WD_HSI_卫图_Level_15.tif"  # Reference image
im_tgt = "D:/LASAC/georeference/ZY1E_AHSI_E29.48_N40.60_20230823_020679_L1A0000644199_SW_EPSG3857.tif"  # Target image to be corrected

# Output paths
output_expanded = r"C:\Users\lyc\OneDrive\Desktop\georeference-output\expanded.tif"  # Expanded image
output_corrected = r"C:\Users\lyc\OneDrive\Desktop\georeference-output\corrected_expanded.tif"  # Corrected image

# Get the bounding box of an image
def get_image_extent(image_path):
    """
    Get the bounding box of an image (coordinates of the bottom-left and top-right corners).
    """
    dataset = gdal.Open(image_path)
    gt = dataset.GetGeoTransform()
    
    # Calculate the coordinates of the four corners of the image
    x_min = gt[0]
    y_max = gt[3]
    x_max = gt[0] + gt[1] * dataset.RasterXSize
    y_min = gt[3] + gt[5] * dataset.RasterYSize
    
    return (x_min, y_min, x_max, y_max)

# Get the bounding boxes of the reference and target images
ref_extent = get_image_extent(im_ref)  # Bounding box of the reference image
tgt_extent = get_image_extent(im_tgt)  # Bounding box of the target image

# Calculate the union of two bounding boxes
def expand_bounding_box(extent1, extent2):
    """
    Calculate the maximum extent of the union of two bounding boxes.
    """
    x_min = min(extent1[0], extent2[0])
    y_min = min(extent1[1], extent2[1])
    x_max = max(extent1[2], extent2[2])
    y_max = max(extent1[3], extent2[3])
    return (x_min, y_min, x_max, y_max)

# Get the maximum extent of the union
expanded_extent = expand_bounding_box(ref_extent, tgt_extent)

# Use gdal.Warp to expand the bounding box of the target image
def expand_image_with_warp(input_path, output_path, output_bounds):
    """
    Use gdal.Warp to expand the bounding box of an image.
    """
    ds = gdal.Open(input_path, gdal.GA_ReadOnly)
    if ds is None:
        raise ValueError("Unable to open the input image!")
    
    # Define Warp options
    warp_options = gdal.WarpOptions(
        format='GTiff',
        outputBounds=output_bounds,  # Expanded bounding box
        dstNodata=0  # Fill value for invalid areas
    )
    
    # Perform the Warp operation
    gdal.Warp(output_path, ds, options=warp_options)
    
    # Close the dataset
    ds = None
    print(f"Expanded image saved to: {output_path}")

# Expand the bounding box of the target image
expand_image_with_warp(im_tgt, output_expanded, expanded_extent)

# Use arosics for automatic correction
coreg = COREG(
    im_ref, 
    output_expanded,  # Use the expanded image as input
    max_shift=1000, 
    ws=(2048, 2048), 
    path_out=output_corrected
)

# Correct the image
coreg.correct_shifts()

print(f"Corrected image saved to: {output_corrected}")

[danschef]

Thanks, it would be helpful to have a tiny code snippet that just reproduces the issue described above. So far as I understand, your posts only contain the STDOUT output of AROSICS and your approach to workaround the issue. If the issue of the missing image areas persists, I would like to reproduce it somehow.

[LiamDiane]

Thank you, I understand what you mean. Here is the code that reproduces the issue mentioned above:

from arosics import COREG
from osgeo import gdal
from shapely.geometry import Polygon

# Input image paths
im_ref = "D:/LASAC/georeference/WD_HSI_卫图/WD_HSI_卫图_Level_15.tif"
im_tgt = "D:/LASAC/georeference/ZY1E_AHSI_E29.48_N40.60_20230823_020679_L1A0000644199_SW_EPSG3857.tif"
# im_tgt = "C:/Users/lyc/OneDrive/Desktop/georeference-output/test.tif"
# Perform image registration
coreg = COREG(
    im_ref, 
    im_tgt, 
    max_shift=1000, 
    ws=(2048, 2048), 
    path_out="C:/Users/lyc/OneDrive/Desktop/georeference-output/test_auto.tif",
)

# Correct the image
coreg.correct_shifts()

print("success")

The AROSICS standard output and related images I shared earlier are the results generated by this code.

[danschef]

Thanks. Would you be able to share the data? The first band from each image would be enough.

[LiamDiane]

How should I send it to you, by email? Is it okay if the compressed file is about 2GB?

[danschef]

Yes, just upload it somewhere and send me the link.

[LiamDiane]

OK.
https://drive.google.com/file/d/1YEknGGCBVWUXOhQHgyIEgmoyeKACQ6zM/view?usp=sharing