Add method to apply coreg to ortho-image

[adehecq]

It's the second time I've been asked about this in less than a week so I raise an issue so we keep this in mind.

Sometimes, people have a DEM and ortho-image of an area. They coregister the DEM to a reference but would like to apply the same transformation to the ortho-image. I don't think we have a straightforward way to do this. I see two options: 1) the easiest solution at the moment is to get the horizontal shift from the coreg.meta dictionary and apply it to the image with Raster.shift. It's simple but it would need to be documented so people can find this information and it only works for a shift (not for more complex affine transformation). 2) in theory, it would be possible to apply the same transformation to the image. However two things need to be taken into account: a) the image may be on a different grid than the DEM, especially different pixel resolution b) any vertical shift should not be applied to the image array directly. Hence one cannot just call Coreg.apply.

Additionally, we would need a test sample of associated DEM and ortho-image for developing, documenting and testing the functionality, which we don't have at the moment.

[rhugonnet]

Good idea. Indeed for 1 the solution would likely be better documentation (now the function is Raster.translate!).

Another option could be to have a short Coreg.apply_horizontal_translation for all AffineCoreg methods, which can be applied to any raster? This way the user wouldn't have to worry about the direction of the shift extracted from the Coreg.meta.

[erikmannerfelt]

I had a similar problem and solution in our terrA project!

My approach was perhaps not the most computationally cheapest, nor maybe the easiest, but here is a solution for a fully 3D approach:

1. Resample the DEM and ortho to the same grid
2. Convert the two products to a colored point cloud (XYZRGB or XYZL where L is lightness)
3. Apply an affine transformation to the point cloud
4. Re-grid the RGB/L part of the point cloud and save.

https://github.com/VAW-SwissTerra/SwissTerra/blob/5b8c2cdf094162c58afe8b488ca1fa19bf919b3f/terra/orthomosaics.py#L29-L39

For simple X/Y shifts, of course something much simpler should be done than this, such as what @rhugonnet just suggested!

I've been thinking of using skimage.transform.warp() for a more general high(-er) performance 3D alternative. It would be something like this:

1. Convert the DEM to a point cloud and add an index to each point.
2. Transform the DEM and record the resultant X/Y shift (in meters) for each point (using the index to match; perhaps the order of the array is enough and would be a faster matching alternative).
3. Convert these X/Y (meters) shifts to X/Y in pixels, and then use skimage.transform.warp() with these shifts to calculate source/destination pixel coordinates of the raster.

I think this would be much easier and could optionally be sped up by subsampling the DEM.

[rhugonnet]

Thanks @erikmannerfelt, I didn't realize we should also consider 3D transformations here, I was thinking only in 2D! Is this what you meant @adehecq? I'm not sure I fully grasp how it makes sense to 3D transform an ortho-image relative to a DEM.

[erikmannerfelt]

Thanks @erikmannerfelt, I didn't realize we should also consider 3D transformations here, I was thinking only in 2D! Is this what you meant @adehecq? I'm not sure I fully grasp how it makes sense to 3D transform an ortho-image relative to a DEM.

If you're asking what the use is, topography makes the transformation non-affine as remapping to different XYZ coordinates is needed.