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))
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.