cmosig / sentle

Sentinel-1 & Sentinel-2 data cubes at large scale (bigger-than-memory) on any machine with integrated cloud detection, snow masking, harmonization, merging, and temporal composites.
MIT License
23 stars 2 forks source link

Possibly unstable georeferencing #26

Closed msoechting closed 1 month ago

msoechting commented 1 month ago

When using Sentle to generate a cube in my test case, the georeferencing seems to be not super stable. This can be seen as a "shaky" animation when putting the cube into Lexcube and animating through time. The area in this case is in the intersection area of two Sentinel tiles, which could be the only reason or just an amplification of an underlying problem.

Cube Generation:

from sentle import sentle
from rasterio.crs import CRS

if __name__ == "__main__":

    da = sentle.process(
            target_crs=CRS.from_string("EPSG:8857"),
            bound_left=961600,
            bound_bottom=6122140,
            bound_right=963070,
            bound_top=6123050,
            datetime="2016-06-17/2024-06-17",
            target_resolution=10,
            S2_mask_snow=True,
            S2_cloud_classification=True,
            S2_cloud_classification_device="cuda",
            S1_assets=None,
            S2_apply_snow_mask=True,
            S2_apply_cloud_mask=True,
            time_composite_freq=None,
            num_workers=30,
        )

    sentle.save_as_zarr(da, "test_cube.zarr")

Visualization:

import lexcube
import numpy as np
import xarray as xr
import spyndex
import cubo
ds = xr.open_dataset("test_cube.zarr", chunks={}, engine="zarr")["sentle"]

d = spyndex.computeIndex(
    index = ["kNDVI", "NDMI"],
    params = {
        "N": ds.sel(band = "B08"),
        "R": ds.sel(band = "B04"),
        "L": 0.5,
        "S1": ds.sel(band = "B11"),
        "kNN" : 1.0,
        "kNR" : spyndex.computeKernel(
            kernel = "RBF",
            params = {
                "a": ds.sel(band = "B08"),
                "b": ds.sel(band = "B04"),
                "sigma": ds.sel(band = ["B08","B04"]).mean("band")
            }),
    }
)

def filter(a):

    nan_percentage = a.isnull().mean(dim=['y', 'x'])
    valid_times = nan_percentage < 0.4
    return a.where(valid_times.compute(), drop=True)

kndvi = filter(d[0])
ndmi = filter(d[1])

w = lexcube.Cube3DWidget(kndvi, vmin=0.0, vmax=1, cmap="tempo")
w

Result:

https://github.com/user-attachments/assets/a72cc67e-d9fd-4ff4-90cc-624375770196

cmosig commented 1 month ago

Looking at the raw sentinel2 geotiffs that we get from planetary computer, we get the exact same spatial movement. Luckily not a bug here :)