danschef
commented
2 years ago Hi @HRodenhizer,
It appears to me that the default of the correct_shifts method is to apply a single correction to all of the bands, although I am not completely sure this is what's happening when I follow the example here.
Correct, by default, AROSICS assumes that there are no remarkable shifts between multiple bands of the same image. Thus, the misregistration is detected in a single band and then corrected in all bands based on the detected tie points.
How would you recommend to loop through the bands of an image? Is there a method for this that I am missing?
Currently, there is no method to apply a co-registration on each band individually. However, it is not a big deal to achieve that by wrapping a bit of code around the actual call of AROSICS. You can use GeoArray to get single bands from a multispectral image (or alternatively pass GDAL VRT files). Then just compute the co-registered result for each band and finally combine the bands into a multispectral image. Maybe this helps:
from geoarray import GeoArray
from arosics import COREG
pr = '/path/to/ref.tif'
pt = '/path/to/tgt.tif'
ref_b3 = GeoArray(pr).get_subset(zslice=slice(2,3)) # get the third band of the reference image
tgt_b3 = GeoArray(pt).get_subset(zslice=slice(2,3)) # get the third band of the target image
# pass the single bands to AROSICS and make sure the output is always resampled to the same pixel grid
# (you may also use COREG_LOCAL here)
CR_b3 = COREG(ref_b3, tgt_b3, target_xyGrid=[[0, 30], [0, 30]])
CR_b3.calculate_spatial_shifts()
output = CR.correct_shifts()
b3_coreg = output['GeoArray_shifted'] # GeoArray instanceYou may call this in a loop for all bands and then either save the co-registered bands to disk and then stack them using GDAL or use numpy + geoarray for this:
import numpy as np
# make sure all of them cover the same area
b1_geoarray_sub = b1_coreg.get_mapPos((xmin, ymin, xmax, ymax)) # pass map bounds here
b2_geoarray_sub = b2_coreg.get_mapPos((xmin, ymin, xmax, ymax)) # pass map bounds here
# stack the bands and save the output image to disk
stacked_array = np.dstack((b1_geoarray_sub, b2_geoarray_sub))
stack = GeoArray(stacked_array, geotransform=b1_coreg.geotransform, projection=b1_coreg.projection)
stack.save('/path/to/stacked_output.tif', fmt='GTiff')
HRodenhizer
Description
In the paper introducing this package, it states: "If co-registration targets a multi-band data cube, the algorithm can be either applied band-by-band in a loop, or the co-registration results of a single band can be applied to the remaining bands."
I would like to loop over the bands of my input image, but I do not see a built in option to do so. It appears to me that the default of the correct_shifts method is to apply a single correction to all of the bands, although I am not completely sure this is what's happening when I follow the example here. Can you confirm or explain what is happening, if that's not correct?
How would you recommend to loop through the bands of an image? Is there a method for this that I am missing?
If not, the most obvious, but messy, solution I see is to create a loop to initiate a COREG object, calculate the shifts, apply the shifts, and save the output file for each band separately. However, this means I would get several output files in which I would only need one of the bands, so I would have to combine the correct bands from the correct files into one image after the initial processing. With one image I don't think that's really a big deal, but with many images, it's definitely not ideal. I think it would also be possible to not save the output file to disk each time, but to copy the shifted GeoArray into a variable in the environment. However, this still has the issue of creating a bunch of GeoArrays in which there is only one band of interest (and I'm not sure yet how to neatly compile the correct bands).
Neither of these options seem great to me, so I'm hoping you have a recommended method that allows the correction to be applied only to the band of interest.