Hi,
I want to perform time series analysis using a harmonic model. Harmonic features and functions are already incorporated into the data cube. My next step involves applying linear regression to all the pixels and extracting the coefficients to fit the curve. However, I couldn’t find a “linear regression” function. I attempted it with process ID but need some suggestions and UDFs to reference.
Apply harmonic function to the data cube
def apply_harmonic_features(ndvi_cube , harmonic_frequencies):
harmonic_cubes = [ndvi_cube ]
for f in harmonic_frequencies:
# Applying cosine
cos_cube = ndvi_cube .apply_dimension(
process=lambda data: cos(data * f),
dimension='t'
).rename_labels(dimension="bands", target=["cos_{}".format(f)])
# Applying sine
sin_cube = ndvi_cube .apply_dimension(
process=lambda data: sin(data * f),
dimension='t'
).rename_labels(dimension="bands", target=["sin_{}".format(f)])
# Combine the cubes
harmonic_cubes.append(cos_cube)
harmonic_cubes.append(sin_cube)
# Merge all harmonic cubes together
from functools import reduce
harmonic_cube = reduce(lambda x, y: x.merge_cubes(y), harmonic_cubes)
return harmonic_cube
harmonic_cube = apply_harmonic_features(ndvi_cube , harmonic_frequencies)
udf = openeo.UDF(
“”"
import numpy as np
import xarray as xr
from openeo.udf import XarrayDataCube
from sklearn.linear_model import LinearRegression
def apply_datacube(cube: XarrayDataCube, context: dict) → XarrayDataCube:
array = cube.get_array()
harmonic_frequencies = context.get(‘harmonic_frequencies’, )
independent_vars = [f"cos_{i}" for i in harmonic_frequencies] + [f"sin_{i}" for i in harmonic_frequencies]
dependent_var = “NDVI”
X = np.vstack([array[ivar].values.flatten() for ivar in independent_vars]).T
y = array[dependent_var].values.flatten()
# linear regression
model = LinearRegression()
model.fit(X, y)
coefficients = model.coef_
intercept = model.intercept_
# dataarrays for coefficients and intercept
coef_array = xr.DataArray(coefficients, dims=["variables"], coords={"variables": independent_vars})
intercept_array = xr.DataArray(intercept, dims=["intercept"], coords={"intercept": ["intercept"]})
#>> results
result_cube = XarrayDataCube({
"coefficients": coef_array,
"intercept": intercept_array
})
return result_cube
“”"
)
ndvi_coeff = harmonic_cube.apply_dimension(code=udf, dimension=“t”)