Georeference and warp a drone image based on sensor orientation and location using transform.from_gcps()


Patrick Gray
 

I have been flying some drone surveys over the ocean and need to properly project and georeference and warp my images. I have all the information I think I need: lat, lon, altitude, yaw, pitch, and roll along with sensor optics params. I have to imagine this can be done with some existing package and I had thought some rasterio functions were suitable, but I can't find any. So I've been using the cameratransform package to get the GPS positions of the corners of my image to feed into rasterio:
 
import cameratransform as ct
 
# camera parameters
 
cam = ct.Camera(ct.RectilinearProjection(focallength_mm=f,
                                         sensor=sensor_size,
                                         image=image_size),
               ct.SpatialOrientation(elevation_m=img.altitude,
                                     tilt_deg=pitch+sensor_offset,
                                     roll_deg=roll,
                                    heading_deg=yaw))
 
# gps pts are lat lon
cam.setGPSpos(img.latitude, img.longitude, img.altitude)
 
# these are the coordinates of the image corners
coords = np.array([cam.gpsFromImage([0 , 0]), \
    cam.gpsFromImage([image_size[0]-1 , 0]), \
    cam.gpsFromImage([image_size[0]-1, image_size[1]-1]), \
    cam.gpsFromImage([0 , image_size[1]-1])])

From there I'm not certain what to do. I've tried basically saying these corners are GCPs and trying to warp the image in rasterio:
 
import rasterio
 
gcp1 = rasterio.control.GroundControlPoint(row=0, col=0, x=coords[0,1], y=coords[0,0], z=coords[0,2], id=None, info=None)
gcp2 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=0, x=coords[1,1], y=coords[1,0], z=coords[1,2], id=None, info=None)
gcp3 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=image_size[1]-1, x=coords[2,1], y=coords[2,0], z=coords[2,2], id=None, info=None)
gcp4 = rasterio.control.GroundControlPoint(row=0, col=image_size[1]-1, x=coords[3,1], y=coords[3,0], z=coords[3,2], id=None, info=None)
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
                    
    # create rasterio transform
    tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
    
    # I also tried this function but to no avail
    #tsfm = rasterio.warp.calculate_default_transform(rasterio.crs.CRS({"init": "epsg:4326"}), rasterio.crs.CRS({"init": "epsg:4326"}), img.size()[0], img.size()[1], gcps=[gcp1,gcp2,gcp3,gcp4])
 
    crs = rasterio.crs.CRS({"init": "epsg:4326"})
 
    profile.update(
        dtype=rasterio.uint8,
        transform = tsfm,
        crs=crs)
        
    with rasterio.open('example.tif', 'w', **profile) as dst:
        dst.write(src.read().astype(rasterio.uint8), 1)

This does produce an image but it is not properly warped. It is just in the correct location. Since my sensor is a 40 deg off nadir the warping should be reasonably significant, when I warp it with cameratransform's cam.getTopViewOfImage() function I get this:

example image

So I'm curious if I'm missing some rasterio function or ability, isn't what the GCPs should be doing in the transform.from_gcps() function?
 


Sean Gillies
 

Hi,

On Mon, Jan 11, 2021 at 11:10 AM <pgrayobx@...> wrote:
I have been flying some drone surveys over the ocean and need to properly project and georeference and warp my images. I have all the information I think I need: lat, lon, altitude, yaw, pitch, and roll along with sensor optics params. I have to imagine this can be done with some existing package and I had thought some rasterio functions were suitable, but I can't find any. So I've been using the cameratransform package to get the GPS positions of the corners of my image to feed into rasterio:
 
import cameratransform as ct
 
# camera parameters
 
cam = ct.Camera(ct.RectilinearProjection(focallength_mm=f,
                                         sensor=sensor_size,
                                         image=image_size),
               ct.SpatialOrientation(elevation_m=img.altitude,
                                     tilt_deg=pitch+sensor_offset,
                                     roll_deg=roll,
                                    heading_deg=yaw))
 
# gps pts are lat lon
cam.setGPSpos(img.latitude, img.longitude, img.altitude)
 
# these are the coordinates of the image corners
coords = np.array([cam.gpsFromImage([0 , 0]), \
    cam.gpsFromImage([image_size[0]-1 , 0]), \
    cam.gpsFromImage([image_size[0]-1, image_size[1]-1]), \
    cam.gpsFromImage([0 , image_size[1]-1])])

From there I'm not certain what to do. I've tried basically saying these corners are GCPs and trying to warp the image in rasterio:
 
import rasterio
 
gcp1 = rasterio.control.GroundControlPoint(row=0, col=0, x=coords[0,1], y=coords[0,0], z=coords[0,2], id=None, info=None)
gcp2 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=0, x=coords[1,1], y=coords[1,0], z=coords[1,2], id=None, info=None)
gcp3 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=image_size[1]-1, x=coords[2,1], y=coords[2,0], z=coords[2,2], id=None, info=None)
gcp4 = rasterio.control.GroundControlPoint(row=0, col=image_size[1]-1, x=coords[3,1], y=coords[3,0], z=coords[3,2], id=None, info=None)
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
                    
    # create rasterio transform
    tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
    
    # I also tried this function but to no avail
    #tsfm = rasterio.warp.calculate_default_transform(rasterio.crs.CRS({"init": "epsg:4326"}), rasterio.crs.CRS({"init": "epsg:4326"}), img.size()[0], img.size()[1], gcps=[gcp1,gcp2,gcp3,gcp4])
 
    crs = rasterio.crs.CRS({"init": "epsg:4326"})
 
    profile.update(
        dtype=rasterio.uint8,
        transform = tsfm,
        crs=crs)
        
    with rasterio.open('example.tif', 'w', **profile) as dst:
        dst.write(src.read().astype(rasterio.uint8), 1)

This does produce an image but it is not properly warped. It is just in the correct location. Since my sensor is a 40 deg off nadir the warping should be reasonably significant, when I warp it with cameratransform's cam.getTopViewOfImage() function I get this:

example image

So I'm curious if I'm missing some rasterio function or ability, isn't what the GCPs should be doing in the transform.from_gcps() function?

First of all, I apologize for not responding sooner.

You must use the rasterio.warp.reproject function to warp. Warping is never implicit in the dataset read/write methods, these are very low-level and always operate strictly in the dataset's row-column space.


--
Sean Gillies


Patrick Gray
 

Hi Sean and thanks for the reply! I've changed my code to use rasterio.warp.reproject, but the final geotiff is still a rectangle rather than a trapezoid and not actually warping to the GCPs. It is roughly where it should be but not even in the middle of the GCPs which are at each corner of the image. Should I be using that rasterio.transform.from_gcps() function to get the transform for the reproject() function? I must be missing something... Here is my code:

import rasterio
from rasterio.warp import reproject
from rasterio.enums import Resampling
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
        
        print(profile)
        # create the transform from the GCPs
        tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
        print(np.array(src.read()))
        
        src_crs = "EPSG:4326"  # This is the crs of the GCPs
 
        # Destination: a 1024 x 1024 dataset in WGS84
        dst_raster = np.zeros((1024, 1024), dtype=np.uint16)
        dst_crs = "EPSG:4326"
 
        tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
 
        new_raster, dst_transform = reproject(
            source=np.array(src.read()),
            destination=dst_raster,
            dst_transform=tsfm,
            gcps=[gcp1,gcp2,gcp3,gcp4],
            src_crs=src_crs,
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
            #**kwargs
        )
 
        profile.update(
            dtype=rasterio.uint16,
            transform = dst_transform,
            crs=dst_crs,
            width=dst_raster.shape[0], # TODO unsure if this is correct order
            height=dst_raster.shape[1]
        )
 
        print(profile)
 
        with rasterio.open('example_uint16_tsfm.tif', 'w', **profile) as dst:
            dst.write(src.read(1).astype(rasterio.uint16), 1)


And here is what the image looks like in WGS84: