Skip to content

Commit 6bc4788

Browse files
committed
updated the RF udp
1 parent 5110101 commit 6bc4788

File tree

4 files changed

+962
-157
lines changed

4 files changed

+962
-157
lines changed
Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
import openeo
2+
from openeo import Connection
3+
from openeo.extra.spectral_indices import compute_indices
4+
from openeo.processes import array_create
5+
6+
7+
8+
def s1_features(connection: Connection, date, aoi, reducer):
9+
10+
"""
11+
Preprocess Sentinel-1 SAR data by applying backscatter correction,
12+
computing VH/VV ratio and log transformations, then reducing over time.
13+
14+
Args:
15+
connection (Connection): An openEO connection.
16+
date (Tuple[str, str]): Temporal extent as (start_date, end_date)
17+
aoi (dict): Spatial extent
18+
reducer (Any): Reducer (first, last, mean, median)
19+
20+
Returns:
21+
DataCube: The processed and temporally reduced Sentinel-1 data cube.
22+
"""
23+
# load S2 pre-collection
24+
s1_cube = connection.load_collection(
25+
"SENTINEL1_GRD",
26+
temporal_extent=date,
27+
spatial_extent=aoi,
28+
bands=["VH", "VV"]
29+
)
30+
31+
# apply SAR backscatter processing to the collection
32+
s1_cube = s1_cube.sar_backscatter(coefficient="sigma0-ellipsoid")
33+
34+
# apply band-wise processing to create a ratio and log-transformed bands
35+
s1_cube = s1_cube.apply_dimension(
36+
dimension="bands",
37+
process=lambda x: array_create(
38+
[
39+
30.0 * x[0] / x[1], # Ratio of VH to VV
40+
30.0 + 10.0 * x[0].log(base=10), # Log-transformed VH
41+
30.0 + 10.0 * x[1].log(base=10), # Log-transformed VV
42+
]
43+
),
44+
)
45+
46+
s1_cube = s1_cube.rename_labels("bands", ["ratio"] + s1_cube.metadata.band_names)
47+
48+
# scale to int16
49+
s1_cube = s1_cube.linear_scale_range(0, 30, 0, 30000)
50+
51+
return s1_cube.reduce_dimension(reducer=reducer, dimension="t")
52+
53+
def s2_features(connection: Connection, date, aoi, reducer):
54+
"""
55+
Preprocess Sentinel-2 data by loading relevant bands, applying scaling,
56+
and reducing over time using a specified reducer.
57+
58+
Args:
59+
connection (Connection): An openEO connection.
60+
date: Temporal extent as (start_date, end_date)
61+
aoi: Spatial extent
62+
reducer (Any): Reducer (first, last, mean, median)
63+
64+
Returns:
65+
DataCube: The processed and temporally reduced Sentinel-2 datacube.
66+
"""
67+
# load S2 pre-collection
68+
s2_cube = connection.load_collection(
69+
"SENTINEL2_L2A",
70+
temporal_extent=date,
71+
spatial_extent=aoi,
72+
bands=["B02", "B03", "B04", "B08","B12"],
73+
max_cloud_cover=80,
74+
)
75+
76+
scl = connection.load_collection(
77+
"SENTINEL2_L2A",
78+
temporal_extent=date,
79+
spatial_extent=aoi,
80+
bands=["SCL"],
81+
max_cloud_cover=80,
82+
)
83+
84+
mask = scl.process(
85+
"to_scl_dilation_mask",
86+
data=scl,
87+
kernel1_size=17,
88+
kernel2_size=77,
89+
mask1_values=[2, 4, 5, 6, 7],
90+
mask2_values=[3, 8, 9, 10, 11],
91+
erosion_kernel_size=3
92+
)
93+
94+
# Create a cloud-free mosaic
95+
masked_cube = s2_cube.mask(mask)
96+
cf_cube = masked_cube.reduce_dimension(reducer=reducer, dimension="t")
97+
98+
# calculate all indices
99+
indices_list = ["NBR", "BAI"]
100+
indices = compute_indices(cf_cube, indices_list)
101+
102+
# calculate texture features
103+
features_udf = openeo.UDF.from_file("features.py")
104+
features = cf_cube.apply_neighborhood(
105+
process=features_udf,
106+
size=[
107+
{"dimension": "x", "value": 128, "unit": "px"},
108+
{"dimension": "y", "value": 128, "unit": "px"},
109+
],
110+
overlap=[
111+
{"dimension": "x", "value": 32, "unit": "px"},
112+
{"dimension": "y", "value": 32, "unit": "px"},
113+
],
114+
)
115+
116+
# combine the original bands with the computed indices,
117+
merged_cube = cf_cube.merge_cubes(indices)
118+
merged_cube = merged_cube.merge_cubes(features)
119+
return merged_cube
Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
import xarray
2+
import numpy as np
3+
from skimage.feature import graycomatrix, graycoprops
4+
from openeo.metadata import CollectionMetadata
5+
6+
7+
def apply_metadata(metadata: CollectionMetadata, context: dict) -> CollectionMetadata:
8+
return metadata.rename_labels(
9+
dimension = "bands",
10+
target = ["contrast","variance","NDFI"]
11+
)
12+
13+
14+
def apply_datacube(cube: xarray.DataArray, context: dict) -> xarray.DataArray:
15+
"""
16+
Applies spatial texture analysis and spectral index computation to a Sentinel-2 data cube.
17+
18+
Computes:
19+
- NDFI (Normalized Difference Fraction Index) from bands B08 and B12
20+
- Texture features (contrast and variance) using Gray-Level Co-occurrence Matrix (GLCM)
21+
22+
Args:
23+
cube (xarray.DataArray): A 3D data cube with dimensions (bands, y, x) containing at least bands B08 and B12.
24+
context (dict): A context dictionary (currently unused, included for API compatibility).
25+
26+
Returns:
27+
xarray.DataArray: A new data cube with dimensions (bands, y, x) containing:
28+
- 'contrast': GLCM contrast
29+
- 'variance': GLCM variance
30+
- 'NDFI': Normalised Difference Fire Index
31+
"""
32+
33+
# Parameters
34+
window_size = 33
35+
pad = window_size // 2
36+
levels = 256 # For 8-bit images
37+
38+
# Load Data
39+
# data = cube.values # shape: (t, bands, y, x)
40+
41+
#first get NDFI
42+
b08 = cube.sel(bands="B08")
43+
b12 = cube.sel(bands="B12")
44+
45+
# Compute mean values
46+
avg_b08 = b08.mean()
47+
avg_b12 = b12.mean()
48+
49+
# Calculate NDFI
50+
ndfi = ((b12 / avg_b12) - (b08 / avg_b08)) / (b08 / avg_b08)
51+
52+
# Padding the image to handle border pixels for GLCM
53+
padded = np.pad(b12, pad_width=pad, mode='reflect')
54+
55+
# Normalize to 0–255 range
56+
img_norm = (padded - padded.min()) / (padded.max() - padded.min())
57+
padded = (img_norm * 255).astype(np.uint8)
58+
59+
# Initialize feature maps
60+
shape = b12.shape
61+
contrast = np.zeros(shape)
62+
variance = np.zeros(shape)
63+
64+
for i in range(pad, pad + shape[0]):
65+
for j in range(pad, pad + shape[1]):
66+
window = padded[i - pad:i + pad + 1, j - pad:j + pad + 1]
67+
68+
# Compute GLCM
69+
glcm = graycomatrix(window, distances=[5], angles=[0], levels=levels, symmetric=True, normed=True)
70+
71+
# Texture features
72+
contrast[i - pad, j - pad] = graycoprops(glcm, 'contrast')[0, 0]
73+
variance[i - pad, j - pad] = np.var(window)
74+
75+
all_texture = np.stack([contrast,variance,ndfi])
76+
# create a data cube with all the calculated properties
77+
textures = xarray.DataArray(
78+
data=all_texture,
79+
dims=["bands", "y", "x"],
80+
coords={"bands": ["contrast","variance","NDFI"], "y": cube.coords["y"], "x": cube.coords["x"]},
81+
)
82+
83+
return textures

algorithm_catalog/vito/random_forest_firemapping/openeo_udp/generate.py

Lines changed: 21 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,12 @@
44
from openeo.api.process import Parameter
55
from openeo.rest.udp import build_process_dict
66

7+
from eo_extractor import s1_features , s2_features
8+
79

810
def generate() -> dict:
911

10-
conn = openeo.connect("openeo.vito.be").authenticate_oidc()
12+
connection = openeo.connect("openeo.vito.be").authenticate_oidc()
1113

1214
# define input parameter
1315
spatial_extent = Parameter.spatial_extent(
@@ -19,21 +21,28 @@ def generate() -> dict:
1921
name = "temporal_extent",
2022
description = "Temporal extent specified as two-element array with start and end date/date-time."
2123
)
22-
# Load s2 bands and set max cloud cover to be less than 10%
23-
s2_bands = conn.load_collection(
24-
collection_id="SENTINEL2_L2A",
25-
spatial_extent=spatial_extent,
26-
temporal_extent=temporal_extent,
27-
bands=["B04", "B08"],
28-
max_cloud_cover=10,
24+
25+
# Load s1 and s2 features
26+
s1_feature_cube = s1_features(connection,temporal_extent,spatial_extent,"median")
27+
s2_feature_cube = s2_features(connection,temporal_extent,spatial_extent,"median")
28+
# Merge the two feature cubes
29+
inference_cube = s2_feature_cube.merge_cubes(s1_feature_cube)
30+
31+
# link to the trained model: this model has an expiry date
32+
model = "https://openeo.vito.be/openeo/1.2/jobs/j-250707130527460a823273791daa4344/results/assets/ZWNjZTlmZWEwNGI4YzljNzZhYzc2YjQ1YjZiYTAwYzIwZjIxMWJkYTQ4NTZjMTRhYTQ0NzViOGU4ZWQ0MzNjZEBlZ2kuZXU=/2c00312531cc55303c14ba68ff2472d4/randomforest.model.tar.gz?expires=1752502120"
33+
34+
# predict of training data
35+
inference = inference_cube.predict_random_forest(
36+
model=model,
37+
dimension="bands"
2938
)
3039

3140
# Build the process dictionary
3241
return build_process_dict(
33-
process_graph=s2_bands,
34-
process_id="s2_bands",
35-
summary="s2_bands",
36-
description="s2_bands",
42+
process_graph=inference,
43+
process_id="random_forest_firemapping",
44+
summary="Forest Fire Mapping Using Random Forest in openEO",
45+
description="Forest fire mapping is a critical tool for environmental monitoring and disaster management, enabling the timely detection and assessment of burned areas. This service is build upon techniques described in the research paper by Zhou, Bao et al., which introduces a machine learning–based approach using Sentinel-2 imagery. Their method combines spectral, topographic, and textural features to improve classification accuracy, particularly emphasising GLCM texture features extracted from Sentinel-2's short-wave infrared band.Thus, the UDP performs forest fire mapping using a pre-trained Random Forest model in openEO. It combines Sentinel-1 and Sentinel-2 features, applies the model, and outputs the predicted fire mapping results.",
3746
parameters=[spatial_extent, temporal_extent]
3847
)
3948

0 commit comments

Comments
 (0)