Image_Lowess.py

import numpy as np
import rasterio as rio
from functools import partial
from statsmodels.nonparametric.smoothers_lowess import lowess


partial_lowess = partial(
    lowess,
    frac=2/3,
    it=3,
    delta=0.0,
    xvals=None,
    is_sorted=False,
    missing='drop',
    return_sorted=True
)

apply_lowess = np.vectorize(partial_lowess, signature="(m),(m)->(m, 2)")

with rio.open("ascat_raw_mask.tif") as ascat_src, rio.open("fac_ais_mask.tif") as fac_src:
    assert ascat_src.crs == fac_src.crs, "CRS mismatch!"
    ascat_data = ascat_src.read(masked=True)
    fac_data = fac_src.read(masked=True)
    ascat_backward = np.array(~ascat_src.transform).reshape((3, 3))
    fac_forward = np.array(fac_src.transform).reshape((3, 3))

    ascat_mask = np.any(ascat_data.mask, axis=0)
    fac_mask = np.any(fac_data.mask, axis=0)

    fac_row, fac_col = np.where(~fac_mask)
    fac_indexes = np.stack(
        [(fac_col + 0.5), (fac_row + 0.5), np.ones_like(fac_col)],
        axis=-1
    )
    ascat_indexes = np.round(
        (
            np.einsum(
                "ij,kj -> ki", ascat_backward, np.einsum(
                    "ij,kj -> ki", fac_forward, fac_indexes
                )
            )[:, :-1] - 0.5
        ),
        decimals=0
    ).astype(int)
    ascat_col, ascat_row = np.array_split(
        ary=ascat_indexes,
        indices_or_sections=ascat_indexes.shape[-1],
        axis=-1
    )

    ascat_col = np.squeeze(a=ascat_col, axis=-1)
    ascat_row = np.squeeze(a=ascat_row, axis=-1)

    indices = ~ascat_mask[ascat_row, ascat_col]
    ascat_col = ascat_col[indices]
    ascat_row = ascat_row[indices]
    fac_row = fac_row[indices]
    fac_col = fac_col[indices]

    ascat_values = ascat_data[:, ascat_row, ascat_col].filled()
    fac_values = fac_data[:, fac_row, fac_col].filled()
    smooth = apply_lowess(ascat_values.T, fac_values.T)