I have Python code that will read drone imagery EXIF data and extract the necessary information (pitch, roll, yaw, elevation...) and leverage that information to create polygons that represent the footprints of those images on the ground. I can manually georeference images to the polygons and confirm that the image footprints are reasonably close to the manually georeferenced images.
I extract the polygon corners and use those coordinates, and the image size, to create GCP arrays that I feed into gdal.gcp to automate the georeferencing the images. The polygon corner order is the same as the image corners starting in the upper left of the image and going clockwise around the image. In the picture below the corners are numbered 0, 1, 2, 3. This is the order that the GCP values are given to GDAL. As you can see from the first image below there is not a perfect alignment of the images to the polygon. I assume there is an affine transformation happening because of the four GCP points. The polygon is not symmetrical.
I have read that GDAL Warp will apply a second order polynomial transfomation (or a thin plate spline) if Warp has 6 GCPs. Hoping to improve the georeferencing to a second order polynomial transformation I have used the polygon corner coordinates to generate the midpoints along the top and bottom of the footprints. Using the image size to get the image midpoint pixels I create two more GCPs and feed those into gdal.gcp. This works without error and the GCPs are now added in order 0,1,2,3,4,5 where 4 and 5 are the top and bottom image edge mid-points.
Here is an example of an images GCPs and image positions in a Python list where each tuple index is the same as each control point. For example, item 0 is corner 0 and the values are (long, lat, imageX, imageY).
[(-122.17137702828244, 48.681497061113696, 0, 0), (-122.16914269018322, 48.682278149382796, 5472, 0), (-122.1684778503216, 48.68129963757258, 5472, 3648), (-122.1705249212658, 48.68058401862344, 0, 3648), (-122.17025985923283, 48.68188760524825, 2736, 0), (-122.1695013857937, 48.680941828098014, 2736, 3648)]
If I add these new georeferenced image outputs to QGIS is see that the georeferencing is somehow corrupt or damaged and the images now display with the upper left image at null island (0, 0). The images are no longer rotated, scaled, or skewed. See picture 2 below.
I have included the code that does the georeferencing below. I can confirm that there is not an error with the two new coordinates, nor the related image coordinates. I do not understand why the addition of these two new points would destroy the georeferencing since I have seen examples online that use more than 4 GCPs with gdal.gcp. I am new to GDAL and the documentation and online examples are few.
How do I apply a second order polynomial transformation using GDAL and my existing code?
Here is the georeferening code.
def rename_image(in_image, out_dir):
'''the output image needs to be a geoTiff so build the path for the file.'''
image_name = in_image.split('/')[-1:]
image_name = image_name[0].split('.')[0] + '.TIF'
out_image = os.path.join(out_dir, image_name)
return out_image
def convert_jpg_to_geotiff(in_image, out_image):
#Convert the jpg to geotiff - this is needed because you need to feed the GCPs to a geoTIFF....
kwargs = {'format': 'GTIFF'}
ds = gdal.Open(in_image)
#Create the geoTiff...
out_file = gdal.Translate(out_image, ds, **kwargs)
#This closes the jpeg file....
ds = None
def build_center_points(coords):
'''creates faux top edge and bottom edge center points to facilitate second order polynomial transformations
This is currently only working for northern hemisphere and western hemisphere becuase of the negative values
for the longitudes and positive values for the lattitudes.'''
tl_long = coords[0][0]
tl_lat = coords[0][1]
tr_long = coords[1][0]
tr_lat = coords[1][1]
bl_long = coords[2][0]
bl_lat = coords[2][1]
br_long = coords[3][0]
br_lat = coords[3][1]
#get geometric centers for top center
tc_lat = (abs(tl_lat) + abs(tr_lat))/2.0
tc_long = -((abs(tl_long) + abs(tr_long))/2.0)
#get geometric centers for bottom center
bc_lat = (abs(bl_lat) + abs(br_lat))/2.0
bc_long = -((abs(bl_long) + abs(br_long))/2.0)
#add ifs to convert to negative values....
coords.append([tc_long, tc_lat])
coords.append([bc_long, bc_lat])
print('top center ', [tc_long, tc_lat])
print('bottom center ',[bc_long, bc_lat])
print(coords)
return coords
def get_center_pixels(image_size):
'''get the top and bottom edge center pixel position'''
image_center = int(image_size[0]/2)
return image_center
def build_gcp_items(coords, image_size, image_center):
'''build the list of list items to feed to the georeferencer'''
#where coords is a list of the image coordinates [(0y,0x), (1y, 1x)...] and image size as a tuple, for example (5472, 3648)
print(coords)
#image index = row, column but gdal calls it GCPPixel, GCPline
gcp_items = [(coords[0][0], coords[0][1], 0, 0), (coords[1][0], coords[1][1], image_size[0], 0), (coords[2][0], coords[2][1], image_size[0], image_size[1]), (coords[3][0], coords[3][1], 0, image_size[1]), (coords[4][0], coords[4][1], image_center, 0), (coords[5][0], coords[5][1], image_center, image_size[1])]
print(gcp_items)
return gcp_items
def georeference(gpc_items, out_image):
'''georeference the geoTiff'''
z = -119
gcp_list = []
counter = 0
for item in gpc_items:
if counter <4:
gcp = gdal.GCP(item[0], item[1], z, item[2], item[3])
print(gcp)
gcp_list.append(gcp)
counter +=1
ds = gdal.Open(out_image, gdal.GA_Update)
wkt = ds.GetProjection()
ds.SetGCPs(gcp_list, wkt)
ds = None
pixel line easting northing.