Topics

Unclear how to implement image alignment based on GCPs

mmweene15@...
 

I am trying to replicated the gdal functionality for performing image alignment based on GCPs. The equivalent gdal commands are the following:

`gdal_translate -of GTiff -gcp GCP1 -gcp GCP2... src_raster dst_raster`
`gdalwarp -s_srs S_SRS -t_srs T_SRS -tps dst_raster dst_raster2`

This would create a new raster, `dst_raster2` that has been reprojected according to the specified GCPs. It is unclear, based on the documentation, how to achieve this using the rasterio.warp module. I have tried using the rasterio.warp.reproject function to no avail. I've included a snippet of code below. The resulting image does not match the image produced by gdal. Any guidance would be appreciated. Note that the variables are descriptively named.

with rasterio.open(target_image_path, 'r') as f:
target_profile = f.profile
target_image = image_tools.rasterio_read_bands(f)
print(target_profile)

destination = np.empty(target_image.shape, dtype=target_profile['dtype'])
resolution = (target_profile['transform'].a, target_profile['transform'].e)
transform, width, height = rasterio_warp.calculate_default_transform(src_crs,
dst_crs,
gcps=gcps,
resolution=resolution,
height=target_profile['height'],
width=target_profile['width'])

rasterio_warp.reproject(source=target_image,
destination=destination,
gcps=gcps,
src_crs=src_crs,
dst_crs=src_crs,
# src_transform=target_profile['transform'],
dst_transform=target_profile['transform'],
resampling=getattr(enums.Resampling, resampling_method),
num_threads=num_threads)

Sean Gillies
 

Hi,

I'm having trouble understanding the intent of code. It looks like you are calling reproject() with source and destination arrays of the same shape, with the same coordinate reference system on each side of the reprojection. I'm not what we should expect in this situation; we are not currently testing for this. Can you write back and show your GCPs?

mmweene15@...
 

Hi, thanks for your response. Let me provide some background. I have two sets of images that capture the same scene but were taken using separate cameras. This causes some degree of GPS discrepancy between the two sets of images, and the task is to align them based on GCPs generated either manually or automatically (irrelevant for this problem). Maybe I have completely misunderstood the purpose of the `reproject` function. I have included a sample of my GCPs below:

`GroundControlPoint(row=array(357), col=array(682), x=array(22.67637871), y=array(-33.94166626), z=0.0, id='45e0b74c-54a4-48c8-b876-1dfcc184c6c8')`

The row and column reference the position in the aligned image that corresponds to the x and y coordinates in the misaligned image. Perhaps my understanding of how to construct this is wrong as well? The documentation doesn't seem to cover this case all that clearly.

Thanks for your help.

Sean Gillies
 

Hi,

On Tue, Nov 6, 2018 at 10:04 PM <mmweene15@...> wrote:
Hi, thanks for your response. Let me provide some background. I have two sets of images that capture the same scene but were taken using separate cameras. This causes some degree of GPS discrepancy between the two sets of images, and the task is to align them based on GCPs generated either manually or automatically (irrelevant for this problem). Maybe I have completely misunderstood the purpose of the `reproject` function. I have included a sample of my GCPs below:

What an interesting problem. The reproject() function is mainly for cartographic reprojection, but it might be useful in your case.
 

`GroundControlPoint(row=array(357), col=array(682), x=array(22.67637871), y=array(-33.94166626), z=0.0, id='45e0b74c-54a4-48c8-b876-1dfcc184c6c8')`

The row and column reference the position in the aligned image that corresponds to the x and y coordinates in the misaligned image. Perhaps my understanding of how to construct this is wrong as well? The documentation doesn't seem to cover this case all that clearly.

The row, col, x, and y for a GroundControlPoint need to be floats, not arrays of floats. I'm sorry that the documentation of this class is inadequate!


Thanks for your help.

You're welcome! I think that adjusting the GCPs may help, but may not be entirely enough.

--
Sean Gillies