Date
1 - 3 of 3
Georeference and warp a drone image based on sensor orientation and location using transform.from_gcps()
|
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:  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:
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
|
|
|
|
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:
|
|
|