How to georeference an unreferenced aerial image using ground control points in python

Question

I have a series of unreferenced aerial images that I would like to georeference using python. The images are identical spatially (they are actually frames extracted from a video), and I obtained ground control points for them by manually georeferencing one frame in ArcMap. I would like to apply the ground control points I obtained to all the subsequent images, and as a result obtain a geo-tiff or a jpeg file with a corresponding world file (.jgw) for each processed image. I know this is possible to do using arcpy, but I do not have access to arcpy, and would really like to use a free open source module if possible.

My coordinate system is NZGD2000 (epsg 2193), and here is the table of control points I wish to apply to my images:

176.412984, -310.977264, 1681255.524654, 6120217.357425

160.386905, -141.487145, 1681158.424227, 6120406.821253

433.204947, -310.547238, 1681556.948690, 6120335.658359

Here is an example image: https://imgur.com/a/9ThHtOz

I've read a lot of information on GDAL and rasterio, but I don't have any experience with them, and am failing to adapt bits of code I found to my particular situation.

Rasterio attempt:

import cv2
from rasterio.warp import reproject
from rasterio.control import GroundControlPoint
from fiona.crs import from_epsg

img = cv2.imread("Example_image.jpg")

# Creating ground control points (not sure if I got the order of variables right):
points = [(GroundControlPoint(176.412984, -310.977264, 1681255.524654, 6120217.357425)),
          (GroundControlPoint(160.386905, -141.487145, 1681158.424227, 6120406.821253)),
          (GroundControlPoint(433.204947, -310.547238, 1681556.948690, 6120335.658359))]

# The function requires a parameter "destination", but I'm not sure what to put there.
#   I'm guessing this may not be the right function to use
reproject(img, destination, src_transform=None, gcps=points, src_crs=from_epsg(2193),
                        src_nodata=None, dst_transform=None, dst_crs=from_epsg(2193), dst_nodata=None,
                        src_alpha=0, dst_alpha=0, init_dest_nodata=True, warp_mem_limit=0)

GDAL attempt:

from osgeo import gdal 
import osr

inputImage = "Example_image.jpg"
outputImage = "image_gdal.jpg"

dataset = gdal.Open(inputImage) 
I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)

outdataset = gdal.GetDriverByName('GTiff') 
output_SRS = osr.SpatialReference() 
output_SRS.ImportFromEPSG(2193) 
outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0]) 
for nb_band in range(I.shape[0]):
    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])

# Creating ground control points (not sure if I got the order of variables right):
gcp_list = [] 
gcp_list.append(gdal.GCP(176.412984, -310.977264, 1681255.524654, 6120217.357425))
gcp_list.append(gdal.GCP(160.386905, -141.487145, 1681158.424227, 6120406.821253))
gcp_list.append(gdal.GCP(433.204947, -310.547238, 1681556.948690, 6120335.658359))

outdataset.SetProjection(srs.ExportToWkt()) 
wkt = outdataset.GetProjection() 
outdataset.SetGCPs(gcp_list,wkt)

outdataset = None

I don't quite know how to make the above code work, and I would really appreciate any help with this.

Answer 1

I ended up reading a book "Geoprocessing with Python" and finally found a solution that worked for me. Here is the code I adapted to my problem:

import shutil
from osgeo import gdal, osr

orig_fn = 'image.tif'
output_fn = 'output.tif'

# Create a copy of the original file and save it as the output filename:
shutil.copy(orig_fn, output_fn)
# Open the output file for writing for writing:
ds = gdal.Open(output_fn, gdal.GA_Update)
# Set spatial reference:
sr = osr.SpatialReference()
sr.ImportFromEPSG(2193) #2193 refers to the NZTM2000, but can use any desired projection

# Enter the GCPs
#   Format: [map x-coordinate(longitude)], [map y-coordinate (latitude)], [elevation],
#   [image column index(x)], [image row index (y)]
gcps = [gdal.GCP(1681255.524654, 6120217.357425, 0, 176.412984, 310.977264),
gdal.GCP(1681158.424227, 6120406.821253, 0, 160.386905, 141.487145),
gdal.GCP(1681556.948690, 6120335.658359, 0, 433.204947, 310.547238)]

# Apply the GCPs to the open output file:
ds.SetGCPs(gcps, sr.ExportToWkt())

# Close the output file in order to be able to work with it in other programs:
ds = None

Answer 2

For your gdal method, just using gdal.Warp with the outdataset should work, e.g.

outdataset.SetProjection(srs.ExportToWkt()) 
wkt = outdataset.GetProjection() 
outdataset.SetGCPs(gcp_list,wkt)
gdal.Warp("output_name.tif", outdataset, dstSRS='EPSG:2193', format='gtiff')

This will create a new file, output_name.tif.