Georeferencing using GDAL geotransform top left corner coordinates

Question

Following this stack exchange example by kflaw, and other examples of the use of geotrasform I am trying to georeference images using pre-computed extent rectangles. However, I am also struggling to correctly position the top left corner. Currently the level of offset increases with capture angle (from north).

Where the capture angle is small, the level of offset is minimal:

What is the correct method for determining the top left coordinates? are they pre or post rotation?

# import sys
# print(sys.path)
from osgeo import gdal, osr
import math

# Camera parameters (constants)
sensorWidth = 36
focalLength = 34.2
imageWidth = 8688
imageHeight = 5792

# Dict of test images
images = {
    # '20180319_091011_029A7639':{ #1
    #   'easting': 1580650.51,
    #   'northing': 5231861.27,
    #   'angle': 105.49,
    #   'altitude': 796
    #   },
    '20180319_092454_029A7649':{ #12
        'easting': 1575143.35,
        'northing': 5224588.46,
        'angle': 346.8,
        'altitude': 715.6,
        'width': 753.31,
        'height': 502.18
    },
    # '20180319_102352_029A7720':{ #83
    #   'easting': 1578520,
    #   'northing': 5212668.253,
    #   'angle': 2.24,
    #   'altitude' : 381.1
    # }
}

# Selected image
image = '20180319_092454_029A7649'

# Set the source and output paths and filenames
src_filename = image+'.tif'
dst_filename = image+'_ref.tif'

# Get the x and y coordinates of the centroid
x = images[image]['easting']
y = images[image]['northing']

# Get the flight vector
inputAngle = images[image]['angle']

# # Compute the GSD as a function of height, image size and camera specifications
# cellSize = ((sensorWidth * images[image]['altitude'] * 100) / (focalLength * imageWidth)) / 100
# print('GSD',cellSize)

cellSizeX =  (images[image]['width']) / imageWidth 
cellSizeY = (images[image]['height']) / imageHeight


# Work out the projected dimensions of the image
width = cellSizeX * imageWidth
height = cellSizeY * imageHeight

# Calculate the x and y scale and skew
angleRadians = math.radians(inputAngle)
scalex = cellSizeX * math.cos(angleRadians)
skewx = -cellSizeX * math.sin(angleRadians)
skewy = -cellSizeY * math.sin(angleRadians)
scaley = -cellSizeY * math.cos(angleRadians)

# Find the coordinates of the upper left corner of the image
upperLeftX = x - width / 2
upperLeftY = y + height / 2
# upperLeftVal = RotateXY(upperLeftX,upperLeftY,x,y,inputAngle)

# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

# Specify raster location through geotransform array
gt = [upperLeftX, scalex, skewx, upperLeftY, skewy, scaley]

# Set location
dst_ds.SetGeoTransform(gt)

# Get raster projection
epsg = 2193
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()

# Set projection
dst_ds.SetProjection(dest_wkt)

# Close the files
dst_ds = None
src_ds = None