Search code examples
pythoncomputer-visiongisgdalgeotiff

Georeference one image using GCPs from another image

I have two GeoTIFFs, one of which is misaligned spatially. I used opencv to do some quick feature matching between the two to automatically identify GCPs existing in both images. How would I go about warping the misaligned image such that it aligns with source/reference layer? I have the geotransforms for each image and the GCPs in both image and geospace.

import cv2
from PIL import Image
import numpy as np
from osgeo import gdal

reference_image = Image.open(r'reference_image.tif')
target_image = Image.open(r'target_image.tif')
reference = np.array(image1)
target = np.array(image2)

sift = cv2.SIFT_create()

keypoints1, descriptors1 = sift.detectAndCompute(reference, None)
keypoints2, descriptors2 = sift.detectAndCompute(target, None)

bf = cv2.BFMatcher()

matches = bf.match(descriptors1, descriptors2)

matches = sorted(matches, key=lambda x: x.distance)

matching_result = cv2.drawMatches(reference, keypoints1, target, keypoints2, matches[:7], None, flags=2)

points1 = [(keypoints1[match.queryIdx].pt[0], (keypoints1[match.queryIdx].pt[1])) for match in matches][0:7]
points2 = [(keypoints2[match.trainIdx].pt[0], (keypoints2[match.trainIdx].pt[1])) for match in matches][0:7]

reference_tif = gdal.Open(r'reference_image.tif')
target_tif = gdal.Open(r'target_image.tif')

reference_gt = reference_tif.GetGeoTransform()
target_gt = target_image.GetGeoTransform()

reference_coords = list()
target_coords = list()

for i in range(len(points1)):
    affine_transform = Affine.from_gdal(*reference_gt)
    x, y = affine_transform * (points1[i][0], points1[i][1])
    transformed_points1.append((x, y))

for i in range(len(points2)):
    affine_transform = Affine.from_gdal(*target_gt)
    x, y = affine_transform * (points2[i][0], points2[i][1])
    transformed_points1.append((x, y))

enter image description here

Solution

  • GeoTIFFs can be georeferenced in two distinct ways: by having GCPs or by having an affine geotransform and an SRS (projection).

    I suppose that you want to set a number of GCPs and then to compute an affine geotransform for a given SRS (projection).

    There is a GDAL utility for doing this and it is called gdalwarp. You can also use it via the Python/C++ API. It is the same utility that is used to reproject images - because this is what you will be doing. The affine geotransform will be computed via a polynomial approximation.

    https://gdal.org/api/python/osgeo.gdal.html