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://i.sstatic.net/roVNi.jpg

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.

Bloomsbury

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

Bloomsbury

Thank you for this answer! I have a similar problem and this helped me understand how I can approach it. I know it is a long shot that you might respond to this, but do you have any idea if it is possible to do this with only one GCP. I also know the resolution of my pixel in meters. Here is my question if you are interested in answering: gis.stackexchange.com/questions/386863/…. – Dimond

@Dimond looks like you have it all figured out, but please let me know if you would still like my input :) – Bloomsbury

@Bloomsbury thanks for the answer, you made my life easier. For anyone that might want to work with jpg instead of tiff this whole procedure works perfectly but you need to take out the gdal.GA_Update option in gdal.Open in order to read jpg images! Cheers – Apc

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.

Willywillynilly

I was initially hanging up on srs.ExportToWkt(), getting the following error: "NameError: name 'srs' is not defined". I just changed this bit of code to output_SRS.ExportToWkt() and now that part seems to run. Now the whole GDAL script runs and produces an "image_gdal.jpg" file, which is not georeferenced and doesn't have any world files associated with it. I added gdal.Warp("output_name.tif", outdataset, dstSRS='EPSG:2193', format='gtiff') as you recommended, but the "output_name.tif" is nowhere to be found - it doesn't seem to be generated. – Bloomsbury

Answer 3

As an addition to @Kat's answer, to avoid quality loss of the original image file and set the nodata-value to 0, the following can be used.

#Load the original file
src_ds = gdal.Open(orig_fn)

#Create tmp dataset saved in memory
driver = gdal.GetDriverByName('MEM')
tmp_ds = driver.CreateCopy('', src_ds, strict=0)

#
# ... setting GCP....
#

# Setting no data for all bands
for i in range(1, tmp_ds.RasterCount + 1):
    f = tmp_ds.GetRasterBand(i).SetNoDataValue(0)

# Saving as file
driver = gdal.GetDriverByName('GTiff')
ds = driver.CreateCopy(output_fn, tmp_ds, strict=0)

Decern