Coordinates of a UAV's image footprint

Question

Good evening, I'm working on a project related to drones and I'm stuck with the following problem:

I need to calculate the actual coordinates (in UTM coordinate system) of each pixel from the drone images. So far what I have managed to do is to calculate the coordinates (in UTM) of each vertex of the image footprint.

Basically, as far as I understand so far the best way to solve the problem is to calculate the transformation matrix from pixels (pixel_x, pixel_y) to that of UTM (latitude_utm, longitude_utm), assuming that the terrain is flat and knowing the following mapping (image the footprint ):

• Top left: (0,0) ----------------------------------> (lat and lng of footprint for Top Left vertex)
• Top Right: (image width, 0) -------------------> (lat and lng of footprint for Top Right vertex)
• Bottom Right: (image width, image height) -> (lat and lng of footprint for Bottom Right vertex)
• Bottom Left: (0,image height) ----------------> (lat and lng of footprint for Bottom Left vertex)

I already tried the below code. But it seems that while it calculates correctly on the upper left peak, it has a large recall on all the others

Note: I only checked against the four peaks because those are the only ones I can evaluate

Input:

# Coordinate traces Manually entered from previous job (Certainly correct)
footprint_coords = [415050.9194898878, 4490804.087897359] # Top Left
                   [415104.8296622897, 4490781.419603585] # Bottom Left
                   [415088.0877967683, 4490885.646303155] # Top Right
                   [415140.5640306050, 4490859.831518373] # Bottom Right
                   
# Read source image
img_src = cv2.imread(IMG_PATH)

# Get source image parameters
img_width, img_height, img_channels = img_src.shape

# Compute image vertex coordinates (in pixels)
src_img_coords = np.array([[0,0]                    # Top Left
                           [0,img_height]           # Bottom Left
                           [img_width,0]            # Top Right
                           [img_width,img_height]]) # Bottom Right

# Get the transformation matrix
project_matrix, _ = cv2.findHomography(src_img_coords, footprint_coords,0)
                                       
# Pre define Array for all pixels coordinates in UTM system
img_utm_coords = np.empty((img_height,img_width), dtype=object)

# Fill the Array
for i in range(img_height):    # rows == height
    for j in range(img_width): # columns == width
        pixel_coords = np.array([j,i,1])
        utm_coords = np.dot(project_matrix, pixel_coords.T)
        img_utm_coords[i,j] = [float(utm_coords[0]),float(utm_coords[1])]
        
# (number of rows, number of columns) == (height, width) == (4000, 3000)
print('UTM array dims:',img_utm_coords.shape)

# Four points
print('TL: ', top_left, ' -> ', img_utm_coords[0,0])
print('TR: ', top_right, ' -> ', img_utm_coords[0,3999])
print('BR: ', bottom_right, '  -> ', img_utm_coords[2999,3999])
print('BL: ', bottom_left, ' -> ', img_utm_coords[2999,0])

Output:

UTM array dims: (3000, 4000)
TL:  [415050.9194898878, 4490804.087897359]  ->  [415050.90624999994, 4490804.0]
TR:  [415088.0877967683, 4490885.646303155]  ->  [415759.75117659935, 4498152.318627857]
BR:  [415140.564030605, 4490859.831518373]   ->  [431890.4374654905, 4672055.155761664]
BL:  [415104.8296622897, 4490781.419603585]  ->  [431181.59253889107, 4664706.837133807]

Answer 1

Using the Rasterio module (based on GDAL), you can use those 4 vertex coordinates in the rasterio.transform.from_gcps() function and you will get the affine transformation matrix for all the pixels in the image.

https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.from_gcps

Then, you can insert that transformation matrix into the tiff file using rasterio:

# Rasterio imports
from rasterio import open as rtopen
from rasterio.transform import from_gcps
from rasterio.control import GroundControlPoint


# Corners control points (for a squared image)
tl = GroundControlPoint(0, 0, left, top)
bl = GroundControlPoint(h_pixel, 0, left, bot)
br = GroundControlPoint(h_pixel, w_pixel, rigth, bot)
tr = GroundControlPoint(0, w_pixel, rigth, top)
gcps = [tl, bl, br, tr]


# Generate the affine transformation matrix
transform = from_gcps(gcps)


# Set the coordinate reference system
epsg = 2062  ## Here you set your own EPSG
crs = 'epsg:{0}'.format(epsg)


# Insert the matrix into the tiff file
with rtopen(tiff_path , 'r+') as tiff:
    tiff.crs = crs
    tiff.transform = transform

I would like to ask you about the calculus for obtaining the vertex coordinates, are they based on camera parameters? Do you have the code to share?