earthdaily / earthdaily/earthdaily-python-client

Register images in time series cube

Open
#12 1 comment 1 reaction 0 assignees View on GitHub
Dominant language
Python
Stars
31
Forks
10
PR merge metrics
No merged PRs in 30d

Description

I've observed pixel shifts of order a couple of pixels on Sentinel 2 and created the following solution which can be integrated/adapted:

```python
from scipy.ndimage import shift as nd_shift
from skimage.registration import phase_cross_correlation
from scipy import fftpack

def low_pass_filter_fourier(image, cutoff_radius):
"""
Apply a low-pass filter to an image in the Fourier domain.

The function performs the Fourier transform on the input image, multiplies
it with a low-pass filter mask, and then performs an inverse Fourier
transform to get the filtered image.

Parameters:
----------
image : np.ndarray
The input 2D image to be filtered.

cutoff_radius : float
The cutoff radius for the low-pass filter in the Fourier domain.
This value should be between 0 and 0.5, where 0.5 corresponds to
the Nyquist frequency.

Returns:
-------
np.ndarray
The filtered image in the spatial domain.

Example:
-------
>>> image = np.random.rand(100, 100)
>>> filtered_image = low_pass_filter_fourier(image, 0.2)
"""
# Fourier Transform
F = fftpack.fftshift(fftpack.fft2(image))

# Generate low-pass filter mask
rows, cols = image.shape
x = np.fft.fftfreq(rows)
y = np.fft.fftfreq(cols)
x, y = np.meshgrid(x, y, indexing='ij')
radius = np.sqrt(x ** 2 + y ** 2)

# Apply low-pass filter
filter_mask = np.fft.fftshift(radius <= cutoff_radius)
F_filtered = F * filter_mask

# Inverse Fourier Transform
image_filtered = np.fft.ifft2(np.fft.ifftshift(F_filtered)).real

return image_filtered

def register_xarrays(xarray1: xr.DataArray, xarray2: xr.DataArray) -> xr.DataArray:
"""Register xarray2 onto xarray1 using the 'red' band for registration."""

array1 = xarray1.sel(band='red').values
array2 = xarray2.sel(band='red').values

array1 = low_pass_filter_fourier(array1, 0.2)
array2 = low_pass_filter_fourier(array2, 0.2)

shift, _, _ = phase_cross_correlation(array1, array2, upsample_factor=100)
print(shift)

corrected_xarray2 = xarray2.copy()

# Apply shift to all bands
for band in xarray2['band'].values:
corrected_xarray2.loc[dict(band=band)] = nd_shift(xarray2.sel(band=band).values, shift[:2])

return corrected_xarray2
```

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.