1

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.

enter image description here

I have read that GDAL Warp will apply a second order polynomial transfomation 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.

enter image description here

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   
3
  • Are the GCPs parameters in right order? GDAL gdal.org/programs/gdal_translate.html wants them as pixel line easting northing. – user30184 May 6 at 22:11
  • GDAL.GCP gdal.org/python/osgeo.gdal.GCP-class.html asks for the parameters to be in a slightly different order than translate. Still - with four GCPs the georeferencing works, it is just the addition of the two non-corner GCPs that seem to cause the trouble. – GBG May 7 at 15:42
  • Would you mind to delete the EXIF info, it is not relevant and makes it harder to read what is essential? – user30184 May 7 at 20:29
1

Both the four GCP images and the six GCP images display with an affine transformation in ArcGIS Pro 2.7, QGIS 3.4, and QGIS 3.10 on Windows 10.

QGIS 3.18 Ubuntu Linux v20.04 is not displaying the six GCP images correctly. The 6 GCP images displayed with the upper left image corner at the equator/greenwich with no rotation, scaling, or resizing. The problem in this case was QGIS 3.18 (Linux) and not the GDAL.GCP tool.

Running gdalwarp -tps infile outfile on the images with 6 GCPs resulted in a new file that displayed the thin plate spline shift.

2
  • Good work. Have you created a QGIS issue? Option -tps performs thin plate splin, not second order polynomial. – user30184 May 7 at 20:38
  • @user30184. Thin plate spline, yes - corrected. I did report my findings to QGIS. – GBG May 7 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.