From e592015759d2b5492c504108c7f9d67228d81662 Mon Sep 17 00:00:00 2001
From: Melanie Weynants <45990429+melwey@users.noreply.github.com>
Date: Mon, 2 Feb 2026 19:33:24 +0100
Subject: [PATCH 1/8] test for multiple indices computation. Not working yet.

---
 Project.toml          |  2 ++
 src/ARCEMEAnalysis.jl | 63 ++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 62 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 6a0b47b..cddc62c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,6 +15,7 @@ Minio = "4281f0d9-7ae0-406e-9172-b7277c1efa20"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
+SpectralIndices = "df0093a1-273d-40bc-819a-796ec3476907"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StringViews = "354b36f9-a18e-4713-926e-db85100087ba"
 YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c"
@@ -33,6 +34,7 @@ Minio = "0.2.2"
 Preferences = "1.5.1"
 ProgressMeter = "1.11.0"
 Proj = "1.9.0"
+SpectralIndices = "0.2.15"
 Statistics = "1.11.1"
 StringViews = "1.3.7"
 YAXArrays = "0.7.0"
diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index 9c85dc1..feaab94 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -9,6 +9,8 @@ import CSV
 using Statistics: mean
 using DataStructures: SortedDict, counter
 using ProgressMeter: @showprogress
+using SpectralIndices: compute_index
+
 include("download.jl")
 
 const arceme_classes = SortedDict(
@@ -29,7 +31,7 @@ const arceme_classes = SortedDict(
 export arceme_cubename, arceme_open, arceme_starttime, arceme_endtime, arceme_eventdate,
     arceme_coordinates, arceme_ndvi, arceme_rgb, arceme_eventlist, arceme_eventpairs, 
     arceme_classes, arceme_landcover, arceme_optical_band_fingerprints, arceme_radar_fingerprints,
-    time_aggregate_fingerprint, arceme_validpairs
+    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral, arceme_test
 
 """
     _arceme_cubenames(;batch="6")
@@ -206,14 +208,69 @@ Compute the NDVI (Normalized Difference Vegetation Index) for the ARCEME data cu
 function arceme_ndvi(ds)
     ndvi = broadcast(ds.B04, ds.B08, ds.cloud_mask, ds.SCL) do b4, b8, cl, scl
         (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
-        fb4 = b4 / typemax(Int16) # not necessary, right?
-        fb8 = b8 / typemax(Int16)
+        fb4 = boa(b4)
+        fb8 = boa(b8)
         (fb8 - fb4) / (fb8 + fb4)
     end
     ds.cubes[:ndvi] = ndvi
     ds
 end
 
+function arceme_test(ds)
+    indices = broadcast(ds.B04, ds.B08, ds.B02, ds.B11, ds.cloud_mask, ds.SCL) do b4, b8, b2, b11, cl, scl
+        (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
+        fb4 = boa(b4)
+        fb8 = boa(b8)
+        S1 = boa(b11); B = boa(b2)
+        ndvi = (fb8 - fb4) / (fb8 + fb4)
+        bri = ((S1 + fb4) - (fb8 + B))/((S1 + fb4) + (fb8 + B))
+        [ndvi, bri]
+    end
+    # @show typeof(indices)
+    # @show size(indices) # size(indices) = (1000, 1000, 146)
+    ds.cubes[:indices] = indices
+    ds
+end
+# does not return a vector it seems
+
+"""
+    arceme_spectral(ds, indices::Vector{String})
+
+Compute the listed indices using SpectralIndices.jl. Not working.
+"""
+function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2") 
+    if platform=="sentinel2" || platform=="sentinel2a" || platform=="sentinel2b"
+        for index in indices
+            tmp = broadcast(ds.cloud_mask, ds.SCL, ds.B01, ds.B02, ds.B03, ds.B04, ds.B08, ds.B05, ds.B06, ds.B07, ds.B11, ds.B12, ds.B09) do cl, scl, b1, b2, b3, b4, b8, b5, b6, b7, b11, b12, b9
+                # BOA
+                A = boa(b1); B = boa(b2); G = boa(b3); R = boa(b4); N = boa(b8)
+                RE1 = boa(b5); RE2 = boa(b6); RE3 = boa(b7)
+                S1 = boa(b11); S2 = boa(b12);  WV = boa(b9)
+                # apply cloud masking here? (or after)
+                (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return repeat([NaN], length(indices))
+                compute_index(index; A, B, G, R, N, RE1, RE2, RE3, S1, S2, WV, L=0.5)
+            end
+            ds.cubes[Symbol(index)] = tmp
+        end
+    elseif platform=="sentinel1"
+        tmp = broadcast(ds.vv, ds.vh) do VV, VH
+            compute_index(indices; VV,VH)
+        end
+    else
+        error("platform $platform is not supported")
+    end
+    ds.cubes[:spectral] = tmp
+    ds
+end
+
+"""
+    boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000)
+
+Compute the radiometric offsets for Sentinel 2 bands to get values at bottom of atmosphere.
+The default values are valid for data processed from baseline 04.00 (January 2022) onwards.
+The entire Sentinel-2 archive in CDSE has been reprocessed and is now available in baseline 05.xx, with consistent offset.
+"""
+boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000) = (band + BOA_ADD_OFFSET) / QUANTIFICATION_VALUE
 
 """
     arceme_rgb(ds)

From a07704f49d6a25a8134be49b9f0d33a3ead6a224 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?M=C3=A9lanie=20Weynants?=
 <45990429+melwey@users.noreply.github.com>
Date: Tue, 3 Feb 2026 16:59:22 +0000
Subject: [PATCH 2/8] delete arceme_test

---
 src/ARCEMEAnalysis.jl | 50 ++++++++++++++++++++++++++-----------------
 1 file changed, 30 insertions(+), 20 deletions(-)

diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index feaab94..426feae 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -31,7 +31,7 @@ const arceme_classes = SortedDict(
 export arceme_cubename, arceme_open, arceme_starttime, arceme_endtime, arceme_eventdate,
     arceme_coordinates, arceme_ndvi, arceme_rgb, arceme_eventlist, arceme_eventpairs, 
     arceme_classes, arceme_landcover, arceme_optical_band_fingerprints, arceme_radar_fingerprints,
-    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral, arceme_test
+    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral
 
 """
     _arceme_cubenames(;batch="6")
@@ -216,23 +216,6 @@ function arceme_ndvi(ds)
     ds
 end
 
-function arceme_test(ds)
-    indices = broadcast(ds.B04, ds.B08, ds.B02, ds.B11, ds.cloud_mask, ds.SCL) do b4, b8, b2, b11, cl, scl
-        (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
-        fb4 = boa(b4)
-        fb8 = boa(b8)
-        S1 = boa(b11); B = boa(b2)
-        ndvi = (fb8 - fb4) / (fb8 + fb4)
-        bri = ((S1 + fb4) - (fb8 + B))/((S1 + fb4) + (fb8 + B))
-        [ndvi, bri]
-    end
-    # @show typeof(indices)
-    # @show size(indices) # size(indices) = (1000, 1000, 146)
-    ds.cubes[:indices] = indices
-    ds
-end
-# does not return a vector it seems
-
 """
     arceme_spectral(ds, indices::Vector{String})
 
@@ -263,6 +246,33 @@ function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2")
     ds
 end
 
+"""
+    arceme_spectral(ds, index::String)
+
+Compute index using SpectralIndices.jl. Working but not when data is actually requested
+"""
+function arceme_spectral(ds, index::String; platform="sentinel2") 
+    if platform=="sentinel2" || platform=="sentinel2a" || platform=="sentinel2b"
+        tmp = broadcast(ds.cloud_mask, ds.SCL, ds.B01, ds.B02, ds.B03, ds.B04, ds.B08, ds.B05, ds.B06, ds.B07, ds.B11, ds.B12, ds.B09) do cl, scl, b1, b2, b3, b4, b8, b5, b6, b7, b11, b12, b9
+            # BOA
+            A = boa(b1); B = boa(b2); G = boa(b3); R = boa(b4); N = boa(b8)
+            RE1 = boa(b5); RE2 = boa(b6); RE3 = boa(b7)
+            S1 = boa(b11); S2 = boa(b12);  WV = boa(b9)
+            # apply cloud masking here? (or after)
+            (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
+            compute_index(index; A, B, G, R, N, RE1, RE2, RE3, S1, S2, WV, L=0.5)
+        end
+        ds.cubes[Symbol(index)] = tmp
+    elseif platform=="sentinel1"
+        tmp = broadcast(ds.vv, ds.vh) do VV, VH
+            compute_index(index; VV,VH)
+        end
+    else
+        error("platform $platform is not supported")
+    end
+    ds
+end
+
 """
     boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000)
 
@@ -279,10 +289,10 @@ Compute the RGB composite for the ARCEME data cube dataset `ds`.
 """
 arceme_rgb(ds) =
     broadcast(ds.B02, ds.B03, ds.B04) do b, g, r
-        m = typemax(Int16)
+        # m = typemax(Int16)
         # RGB(r / m * 4, g / m * 4, b / m * 4)
         # RGB(min(1.0,r / m *4) , min(1.0,g / m *4) , min(1.0,b / m *4))
-        RGB(r / m , g / m , b / m )
+        RGB(boa(r), boa(g), boa(b))
 end
 
 """

From c4b7b822a87f286cc3364709f4f5060eed937eda Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?M=C3=A9lanie=20Weynants?=
 <45990429+melwey@users.noreply.github.com>
Date: Wed, 4 Feb 2026 12:49:52 +0000
Subject: [PATCH 3/8] set httpstore

---
 src/ARCEMEAnalysis.jl | 20 +++++++++++---------
 src/download.jl       | 24 ++++++++++++++++++++----
 2 files changed, 31 insertions(+), 13 deletions(-)

diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index 426feae..f342bcf 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -40,7 +40,7 @@ List all available data cube names in the specified batch stored in the ARCEME S
 """
 function _arceme_cubenames(; batch="ARCEME-DC-6")
     cubenames = if local_cubepath === nothing
-        store = S3Store("$(batch)/", MinioConfig("https://s3.waw3-2.cloudferro.com/swift/v1"))
+        store = S3Store("$(batch)/", MinioConfig(httpstore))
 
         resp = Zarr.cloud_list_objects(store, batch)
         map(split(String(resp), "\n")) do p
@@ -115,22 +115,23 @@ arceme_eventpairs() =
     Iterators.partition(sort(arceme_eventlist(),by=i->(i.dhp_label,i.source)),2) |> collect
 
 """
-    arceme_validpairs(;batch="ARCEME-DC-6", store="https://s3.waw3-2.cloudferro.com/swift/v1")
+    arceme_validpairs(;batch="ARCEME-DC-6")
 
-Get valid pairs of ARCEME events from the data store. Currently only working over HTTP(S)
+Get valid pairs of ARCEME events from the local path (if set with arceme_set_localpath) or the http data store.
+Default httpstore is "https://s3.waw3-2.cloudferro.com/swift/v1". It can be reset with arceme_set_httpstore.
 """
-function arceme_validpairs(;batch="ARCEME-DC-6", store="https://s3.waw3-2.cloudferro.com/swift/v1")
+function arceme_validpairs(;batch="ARCEME-DC-6")
     allpairs = arceme_eventpairs()
     validpairs = if local_cubepath === nothing
-        map(x -> all([Zarr.is_zgroup(Zarr.HTTPStore("$store/$batch/$(arceme_cubename(i))"), "") for i in x]), allpairs)
+        map(x -> all([Zarr.is_zgroup(Zarr.HTTPStore("$httpstore/$batch/$(arceme_cubename(i))"), "") for i in x]), allpairs)
     else
         map(x -> all([isfile(joinpath(local_cubepath, batch, string(arceme_cubename(i), ".zip"))) for i in x]), allpairs)
     end
     allpairs[validpairs]
 end
 
-function arceme_open(event::Event)
-    arceme_open(arceme_cubename(event))
+function arceme_open(event::Event; batch="ARCEME-DC-6")
+    arceme_open(arceme_cubename(event); batch="ARCEME-DC-6")
 end
 
 """
@@ -151,11 +152,12 @@ arceme_landcover(ev::Event) = arceme_landcover(arceme_open(ev))
 """ 
     arceme_open(cubename; batch="ARCEME-DC-6")
 
-Open the specified ARCEME data cube from the S3 bucket.
+Open the specified ARCEME data cube from the local path (if with arceme_set_localpath) or 
+    the httpstore (default "https://s3.waw3-2.cloudferro.com/swift/v1", reset with arceme_set_httpstore).
 """
 function arceme_open(cubename; batch="ARCEME-DC-6")
     if local_cubepath === nothing
-        open_dataset("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/$cubename", force_datetime=true)
+        open_dataset("$httpstore/$batch/$cubename", force_datetime=true)
     else
         open_dataset(joinpath(local_cubepath, batch, string(cubename, ".zip")))
     end
diff --git a/src/download.jl b/src/download.jl
index 2008df3..5770640 100644
--- a/src/download.jl
+++ b/src/download.jl
@@ -3,7 +3,22 @@ using ProgressMeter
 using HTTP, JSON
 using ZipArchives
 using Preferences: @set_preferences!, @load_preference
-export arceme_set_localpath
+export arceme_set_localpath, arceme_set_httpstore
+
+"""
+arceme_set_httpstore(httpstore)
+
+Sets the https store from which to retrieve the cubes.
+If no local preference is set, the default location is used: https://s3.waw3-2.cloudferro.com/swift/v1
+
+Setting this preference needs recompilation.
+"""
+function arceme_set_httpstore(store)
+    @set_preferences!("arceme_httpstore" => store)
+    @warn "Preferences changed! Restart Julia for this change to take effect."
+end
+
+const httpstore = @load_preference("arceme_httpstore", "https://s3.waw3-2.cloudferro.com/swift/v1")
 
 """
     arceme_set_localpath(path)
@@ -12,6 +27,7 @@ Sets the path to a local directory where downloaded copies of the arceme cubes a
 """
 function arceme_set_localpath(path)
     @set_preferences!("arceme_localpath" => path)
+    @warn "Preferences changed! Restart Julia for this change to take effect."
 end
 
 const local_cubepath = @load_preference("arceme_localpath")
@@ -20,17 +36,17 @@ function arceme_download_batch(batch="ARCEME-DC-6")
     if local_cubepath === nothing
         error("You need to set the local cube path first. Please run `arceme_localpath(path)` first.")
     end
-    aresp = HTTP.get("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/",query=Dict("format"=>"json","delimiter"=>"/"))
+    aresp = HTTP.get("$httpstore/$batch/",query=Dict("format"=>"json","delimiter"=>"/"))
     allarrays = map(i->strip(i["subdir"],'/'),JSON.parse(aresp.body))
 
     @showprogress for current_cube in allarrays
-        lresp = HTTP.get("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/",query=Dict("prefix"=>current_cube))
+        lresp = HTTP.get("$httpstore/$batch/",query=Dict("prefix"=>current_cube))
         files_in_cube = split(StringView(lresp.body),"\n")
 
         outfilename = joinpath(local_cubepath,batch,string(current_cube,".zip"))
         ZipWriter(outfilename) do w
             for f in files_in_cube
-                resp = HTTP.get("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/$f");
+                resp = HTTP.get("$httpstore/$batch/$f");
                 f2 = joinpath(splitpath(f)[2:end]...)
                 zip_newfile(w, f2)
                 write(w,resp.body)

From d86956ca2f488ac6ae61ffdba7f192cf40374685 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?M=C3=A9lanie=20Weynants?=
 <45990429+melwey@users.noreply.github.com>
Date: Wed, 4 Feb 2026 12:59:31 +0000
Subject: [PATCH 4/8] fix batch kwarg

---
 src/ARCEMEAnalysis.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index f342bcf..4b8b5e7 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -131,7 +131,7 @@ function arceme_validpairs(;batch="ARCEME-DC-6")
 end
 
 function arceme_open(event::Event; batch="ARCEME-DC-6")
-    arceme_open(arceme_cubename(event); batch="ARCEME-DC-6")
+    arceme_open(arceme_cubename(event); batch=batch)
 end
 
 """
@@ -146,7 +146,7 @@ function arceme_landcover(ds)
         (key=k, class=v, count=count, fraction=count/1000000)
     end
 end
-arceme_landcover(ev::Event) = arceme_landcover(arceme_open(ev))
+arceme_landcover(ev::Event; batch="ARCEME-DC-6") = arceme_landcover(arceme_open(ev, batch=batch))
 
 
 """ 

From 1191a92f8d8ac6f13990c4dec02d4afc707512d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?M=C3=A9lanie=20Weynants?=
 <45990429+melwey@users.noreply.github.com>
Date: Wed, 4 Feb 2026 17:17:28 +0000
Subject: [PATCH 5/8] arceme_spectral error details

---
 src/ARCEMEAnalysis.jl | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index 4b8b5e7..6c1dfd8 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -251,7 +251,15 @@ end
 """
     arceme_spectral(ds, index::String)
 
-Compute index using SpectralIndices.jl. Working but not when data is actually requested
+Compute index using SpectralIndices.jl, e.g. "NDVI". Not Working (when actually requesting the data).
+
+Example:
+validpairs = arceme_validpairs()
+ds_d,ds_dhp = arceme_open.(validpairs[80])
+arceme_spectral(ds_d, "NDVI")
+@time ds_d.NDVI[x=1,y=1,].data[:]
+ERROR: MethodError: no method matching (::XFunction{ARCEMEAnalysis.var"#36#37"{String}, XOutput{Tuple{}, Tuple{}, Int64}, Tuple{}})
+The function `XFunction{ARCEMEAnalysis.var"#36#37"{String}, XOutput{Tuple{}, Tuple{}, Int64}, Tuple{}}(ARCEMEAnalysis.var"#36#37"{String}("NDVI"), XOutput{Tuple{}, Tuple{}, Int64}((), (), 1, Dict{Any, Any}()), (), false)` exists, but no method is defined for this combination of argument types.
 """
 function arceme_spectral(ds, index::String; platform="sentinel2") 
     if platform=="sentinel2" || platform=="sentinel2a" || platform=="sentinel2b"

From 9cb642eae9dde43184eeb9d73a48492e2d1c3c33 Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Tue, 10 Feb 2026 09:57:17 +0000
Subject: [PATCH 6/8] update spectral function

---
 Project.toml            |  4 ++--
 src/ARCEMEAnalysis.jl   | 31 ++++++++++++++++---------------
 src/spectral_helpers.jl | 38 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 56 insertions(+), 17 deletions(-)
 create mode 100644 src/spectral_helpers.jl

diff --git a/Project.toml b/Project.toml
index cddc62c..ae6fe7d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,7 @@ ZipArchives = "49080126-0e18-4c2a-b176-c102e4b3760c"
 [compat]
 CSV = "0.10.15"
 Colors = "0.13.1"
-DataStructures = "0.18.22"
+DataStructures = "0.18.22, 0.19"
 Dates = "1.11.0"
 DimensionalData = "0.29.24"
 HTTP = "1.10.19"
@@ -37,7 +37,7 @@ Proj = "1.9.0"
 SpectralIndices = "0.2.15"
 Statistics = "1.11.1"
 StringViews = "1.3.7"
-YAXArrays = "0.7.0"
+YAXArrays = "0.7.2"
 Zarr = "0.9.5"
 ZipArchives = "2.6.0"
 julia = "1.10"
diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index 6c1dfd8..68b400a 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -9,7 +9,7 @@ import CSV
 using Statistics: mean
 using DataStructures: SortedDict, counter
 using ProgressMeter: @showprogress
-using SpectralIndices: compute_index
+using SpectralIndices: compute_index, SpectralIndices
 
 include("download.jl")
 
@@ -218,25 +218,23 @@ function arceme_ndvi(ds)
     ds
 end
 
+#Helper functions to compute the indices from named tuples for type stability
+compute_indexx(index::SpectralIndices.SpectralIndex{<:Any,B}, params::NamedTuple) where B = index.compute(Float64, params[B]...)
+function listofindices(indices, values)
+    map(indices) do index
+        compute_indexx(index, values)
+    end
+end
+
 """
     arceme_spectral(ds, indices::Vector{String})
 
 Compute the listed indices using SpectralIndices.jl. Not working.
 """
 function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2") 
-    if platform=="sentinel2" || platform=="sentinel2a" || platform=="sentinel2b"
-        for index in indices
-            tmp = broadcast(ds.cloud_mask, ds.SCL, ds.B01, ds.B02, ds.B03, ds.B04, ds.B08, ds.B05, ds.B06, ds.B07, ds.B11, ds.B12, ds.B09) do cl, scl, b1, b2, b3, b4, b8, b5, b6, b7, b11, b12, b9
-                # BOA
-                A = boa(b1); B = boa(b2); G = boa(b3); R = boa(b4); N = boa(b8)
-                RE1 = boa(b5); RE2 = boa(b6); RE3 = boa(b7)
-                S1 = boa(b11); S2 = boa(b12);  WV = boa(b9)
-                # apply cloud masking here? (or after)
-                (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return repeat([NaN], length(indices))
-                compute_index(index; A, B, G, R, N, RE1, RE2, RE3, S1, S2, WV, L=0.5)
-            end
-            ds.cubes[Symbol(index)] = tmp
-        end
+    tmp = if platform == "sentinel2" || platform == "sentinel2a" || platform == "sentinel2b"
+        pl = platform == "sentinel2" ? "sentinel2a" : platform
+        _compute_indices(ds, indices, pl)
     elseif platform=="sentinel1"
         tmp = broadcast(ds.vv, ds.vh) do VV, VH
             compute_index(indices; VV,VH)
@@ -244,7 +242,9 @@ function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2")
     else
         error("platform $platform is not supported")
     end
-    ds.cubes[:spectral] = tmp
+    foreach(pairs(tmp)) do (k,v)
+        ds.cubes[k] = v
+    end
     ds
 end
 
@@ -416,4 +416,5 @@ end
 
 
 include("spatialdebias.jl")
+include("spectral_helpers.jl")
 end #module
\ No newline at end of file
diff --git a/src/spectral_helpers.jl b/src/spectral_helpers.jl
new file mode 100644
index 0000000..fe072f0
--- /dev/null
+++ b/src/spectral_helpers.jl
@@ -0,0 +1,38 @@
+import SpectralIndices as SI
+normalize_s2name(n) = match(r"B\d",n) === nothing ? Symbol(n) : Symbol("B0$(last(n))")
+struct NTWrapper{F,names,C,I} <: Function
+    f::F
+    names::Val{names}
+    consts::C
+    indices::I
+end
+function (ntw::NTWrapper{F,names})(cl,scl,args...) where {F,names} 
+    ntw.f(cl,scl,NamedTuple{names}(args),ntw.consts,ntw.indices) 
+end
+NTWrapper(f,names::Tuple,consts::NamedTuple,indices) = NTWrapper(f,Val(names),consts,indices)
+function inner_compute_indices(cl,scl,bands,consts,indices_tuple)
+    (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return map(_ -> NaN, indices_tuple)
+    bandparams = map(ARCEMEAnalysis.boa,bands)
+    allparams = (; bandparams..., consts...)
+    listofindices(indices_tuple, allparams)
+end
+
+function _compute_indices(ds,indices,platform)
+    indices_tuple = ((SI.indices[k] for k in indices)...,)
+
+    allvars = mapreduce(r->Set(string.(SI._band_names(r))),union!,indices_tuple)
+
+    needed_bands = intersect(allvars,keys(SpectralIndices.bands))
+    needed_constants = intersect(allvars,keys(SpectralIndices.constants))
+    undefined_constants = setdiff(allvars,needed_bands, needed_constants)
+    default_values = (;(Symbol(c)=>SI.constants[c].default for c in needed_constants)...)
+    bands_to_load = (sort([normalize_s2name(SI.bands[b].platforms[platform].band) for b in needed_bands])...,)
+
+    wrapped_function = NTWrapper(inner_compute_indices, (Symbol.(needed_bands)...,),default_values,indices_tuple)
+
+    bandargs = map(b->ds.cubes[b],bands_to_load)
+    output = map(_ -> XOutput(; outtype=Float32), indices_tuple)
+
+
+    return NamedTuple{(Symbol.(indices)...,)}(xmap(wrapped_function, ds.cloud_mask, ds.SCL, bandargs...; output, inplace=false))
+end
\ No newline at end of file

From 3bc6aeb07f3deb115df4f5bc20784613f080375c Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Fri, 20 Feb 2026 14:23:35 +0000
Subject: [PATCH 7/8] Index computation

---
 Project.toml            |   4 ++
 src/ARCEMEAnalysis.jl   | 113 ++++++++++++++++++++++------------------
 src/s1_helpers.jl       |  45 ++++++++++++++++
 src/spatialdebias.jl    |  30 +++++------
 src/spectral_helpers.jl | 104 ++++++++++++++++++++++++++++++++----
 5 files changed, 218 insertions(+), 78 deletions(-)
 create mode 100644 src/s1_helpers.jl

diff --git a/Project.toml b/Project.toml
index ae6fe7d..ec1773e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,9 +9,11 @@ Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
+DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Minio = "4281f0d9-7ae0-406e-9172-b7277c1efa20"
+OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
@@ -28,9 +30,11 @@ Colors = "0.13.1"
 DataStructures = "0.18.22, 0.19"
 Dates = "1.11.0"
 DimensionalData = "0.29.24"
+DiskArrayEngine = "0.3.0"
 HTTP = "1.10.19"
 JSON = "1.4.0"
 Minio = "0.2.2"
+OnlineStats = "1.7.3"
 Preferences = "1.5.1"
 ProgressMeter = "1.11.0"
 Proj = "1.9.0"
diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index 68b400a..b942e9b 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -1,7 +1,7 @@
 module ARCEMEAnalysis
 using Zarr: S3Store, Zarr
 using Minio: MinioConfig
-using YAXArrays: open_dataset, ⊘, xmap, XOutput, YAXArray, YAXArrays, Dataset
+using YAXArrays: open_dataset, ⊘, xmap, XOutput, YAXArray, YAXArrays, Dataset, setchunks, compute_to_zarr, savedataset
 import DimensionalData as DD
 using Colors: RGB
 using Dates: DateTime, Year, Date
@@ -10,8 +10,20 @@ using Statistics: mean
 using DataStructures: SortedDict, counter
 using ProgressMeter: @showprogress
 using SpectralIndices: compute_index, SpectralIndices
+import SpectralIndices as SI
+
+function __init__()
+    #Extend SpectralIndices band definitions with S1 bands
+    s1vv = Dict{String,SI.PlatformBand}("sentinel1" => SI.PlatformBand("sentinel1", "VV", "Vertical-Vertical", 5.55e7, 1e5))
+    s1vh = Dict{String,SI.PlatformBand}("sentinel1" => SI.PlatformBand("sentinel1", "VH", "Vertical-Horizontal", 5.55e7, 1e5))
+
+    SI.bands["VV"] = SI.Band("VV", "Vertical-Vertical", "vv", 5.54e7, 5.56e7, s1vv)
+    SI.bands["VH"] = SI.Band("VH", "Vertical-Horizontal", "vh", 5.54e7, 5.56e7, s1vh)
+end
+
 
 include("download.jl")
+include("s1_helpers.jl")
 
 const arceme_classes = SortedDict(
   0   => "No data",
@@ -27,11 +39,12 @@ const arceme_classes = SortedDict(
   95  => "Mangroves",
   100 => "Moss and lichen",
 )
+lckeymap(k) = ifelse(k > 90, (k + 10), k) ÷ 10 + 1
 
 export arceme_cubename, arceme_open, arceme_starttime, arceme_endtime, arceme_eventdate,
     arceme_coordinates, arceme_ndvi, arceme_rgb, arceme_eventlist, arceme_eventpairs, 
     arceme_classes, arceme_landcover, arceme_optical_band_fingerprints, arceme_radar_fingerprints,
-    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral
+    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral, arceme_kndvi, arceme_radar_db
 
 """
     _arceme_cubenames(;batch="6")
@@ -159,7 +172,16 @@ function arceme_open(cubename; batch="ARCEME-DC-6")
     if local_cubepath === nothing
         open_dataset("$httpstore/$batch/$cubename", force_datetime=true)
     else
-        open_dataset(joinpath(local_cubepath, batch, string(cubename, ".zip")))
+        main_ds = open_dataset(joinpath(local_cubepath, batch, string(cubename, ".zip")))
+        if isfile(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
+            index_ds = open_dataset(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
+            for (k, v) in (index_ds.cubes)
+                main_ds.cubes[k] = v
+            end
+            return main_ds
+        else
+            return main_ds
+        end
     end
 end
 
@@ -209,7 +231,7 @@ Compute the NDVI (Normalized Difference Vegetation Index) for the ARCEME data cu
 """
 function arceme_ndvi(ds)
     ndvi = broadcast(ds.B04, ds.B08, ds.cloud_mask, ds.SCL) do b4, b8, cl, scl
-        (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
+        _is_cloud(cl, scl) && return NaN
         fb4 = boa(b4)
         fb8 = boa(b8)
         (fb8 - fb4) / (fb8 + fb4)
@@ -218,14 +240,6 @@ function arceme_ndvi(ds)
     ds
 end
 
-#Helper functions to compute the indices from named tuples for type stability
-compute_indexx(index::SpectralIndices.SpectralIndex{<:Any,B}, params::NamedTuple) where B = index.compute(Float64, params[B]...)
-function listofindices(indices, values)
-    map(indices) do index
-        compute_indexx(index, values)
-    end
-end
-
 """
     arceme_spectral(ds, indices::Vector{String})
 
@@ -236,9 +250,7 @@ function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2")
         pl = platform == "sentinel2" ? "sentinel2a" : platform
         _compute_indices(ds, indices, pl)
     elseif platform=="sentinel1"
-        tmp = broadcast(ds.vv, ds.vh) do VV, VH
-            compute_index(indices; VV,VH)
-        end
+        _compute_indices(ds, indices, platform)
     else
         error("platform $platform is not supported")
     end
@@ -248,40 +260,9 @@ function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2")
     ds
 end
 
-"""
-    arceme_spectral(ds, index::String)
-
-Compute index using SpectralIndices.jl, e.g. "NDVI". Not Working (when actually requesting the data).
-
-Example:
-validpairs = arceme_validpairs()
-ds_d,ds_dhp = arceme_open.(validpairs[80])
-arceme_spectral(ds_d, "NDVI")
-@time ds_d.NDVI[x=1,y=1,].data[:]
-ERROR: MethodError: no method matching (::XFunction{ARCEMEAnalysis.var"#36#37"{String}, XOutput{Tuple{}, Tuple{}, Int64}, Tuple{}})
-The function `XFunction{ARCEMEAnalysis.var"#36#37"{String}, XOutput{Tuple{}, Tuple{}, Int64}, Tuple{}}(ARCEMEAnalysis.var"#36#37"{String}("NDVI"), XOutput{Tuple{}, Tuple{}, Int64}((), (), 1, Dict{Any, Any}()), (), false)` exists, but no method is defined for this combination of argument types.
-"""
-function arceme_spectral(ds, index::String; platform="sentinel2") 
-    if platform=="sentinel2" || platform=="sentinel2a" || platform=="sentinel2b"
-        tmp = broadcast(ds.cloud_mask, ds.SCL, ds.B01, ds.B02, ds.B03, ds.B04, ds.B08, ds.B05, ds.B06, ds.B07, ds.B11, ds.B12, ds.B09) do cl, scl, b1, b2, b3, b4, b8, b5, b6, b7, b11, b12, b9
-            # BOA
-            A = boa(b1); B = boa(b2); G = boa(b3); R = boa(b4); N = boa(b8)
-            RE1 = boa(b5); RE2 = boa(b6); RE3 = boa(b7)
-            S1 = boa(b11); S2 = boa(b12);  WV = boa(b9)
-            # apply cloud masking here? (or after)
-            (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
-            compute_index(index; A, B, G, R, N, RE1, RE2, RE3, S1, S2, WV, L=0.5)
-        end
-        ds.cubes[Symbol(index)] = tmp
-    elseif platform=="sentinel1"
-        tmp = broadcast(ds.vv, ds.vh) do VV, VH
-            compute_index(index; VV,VH)
-        end
-    else
-        error("platform $platform is not supported")
-    end
-    ds
-end
+_is_cloud(cl, scl) = (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11)))
+
+
 
 """
     boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000)
@@ -302,7 +283,7 @@ arceme_rgb(ds) =
         # m = typemax(Int16)
         # RGB(r / m * 4, g / m * 4, b / m * 4)
         # RGB(min(1.0,r / m *4) , min(1.0,g / m *4) , min(1.0,b / m *4))
-        RGB(boa(r), boa(g), boa(b))
+        RGB(clamp(boa(r), 0, 1), clamp(boa(g), 0, 1), clamp(boa(b), 0, 1))
 end
 
 """
@@ -413,6 +394,38 @@ function time_aggregate_fingerprint(allbands, eventdate, banddim, timeaxis)
     YAXArray((allbands.lc, DD.Ti((fingerprint_timesteps .- 0.5) ./ step_per_year .* 12), banddim), cat(res..., dims=2))
 end
 
+"""
+For every valid event pair creates a and stores data cubes of a list of precomputed vegetation indices. For kNDVI 
+a shared sigma parameter per land cover class is computed.
+"""
+function arceme_create_indexcubes(; indices_s1=["DpRVIVV"], indices_s2=["NDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"])
+
+    for ev in arceme_validpairs()
+
+        ds_pair = arceme_open.(ev)
+
+        foreach(ds_pair) do ds
+            arceme_spectral(ds, indices_s1, platform="sentinel1")
+            arceme_spectral(ds, indices_s2, platform="sentinel2")
+            arceme_radar_db(ds)
+        end
+        ARCEMEAnalysis.arceme_kndvi_pair(ds_pair...)
+
+        fields_to_save = [indices_s1; indices_s2]
+
+        foreach(ds_pair, ev) do ds, event
+            output_base = "$local_cubepath/ARCEME-DC-6-INDICES"
+            name = arceme_cubename(event)
+            indexcube = setchunks(ds[fields_to_save], (500, 500, 25))
+            compute_to_zarr(indexcube, joinpath(output_base, name), overwrite=true)
+            cube2 = setchunks(ds[["vv_db", "vh_db", "kNDVI"]], (500, 500, 25))
+            savedataset(cube2, path=joinpath(output_base, name), append=true)
+            run(Cmd(`zip -0 -r ../$(name).zip .`, dir=joinpath(output_base, name)))
+            rm(joinpath(output_base, name), recursive=true)
+        end
+    end
+end
+
 
 
 include("spatialdebias.jl")
diff --git a/src/s1_helpers.jl b/src/s1_helpers.jl
new file mode 100644
index 0000000..b453071
--- /dev/null
+++ b/src/s1_helpers.jl
@@ -0,0 +1,45 @@
+import Dates: Millisecond
+import YAXArrays: YAXArray
+"""
+    arceme_s1_position(ds)
+
+Analyses the sentinel 1 time stamps of the datasets and groups them into series of time stamps whose difference
+is always a multiple of 6 days, so we can assume that all images with the same group tag are retrieved from the 
+same position. Store the groups in the dataset with the name `s1_postion`.
+"""
+function arceme_s1_position(ds)
+    ts = ds.time_sentinel_1_rtc.val
+    sixdays = 6*24*60*60*1000
+    group_offsets = [0]
+    groups = [[1]]
+    for its in 2:length(ts)
+        groupfound = false
+        for igroup in eachindex(groups)
+            fac = (Millisecond(ts[its]-ts[1]).value-group_offsets[igroup])/sixdays
+            if abs(round(fac)-fac) < 1e-6
+                push!(groups[igroup],its)
+                groupfound=true
+                break
+            end
+        end
+        if !groupfound
+            newoffset = Dates.Millisecond(ts[its]-ts[1]).value-floor(Int,fac)*sixdays
+            push!(group_offsets,newoffset)
+            push!(groups,[its])
+        end
+    end
+    positions = zeros(Int,length(ts))
+    for i in 1:length(groups)
+        positions[groups[i]] .= i
+    end
+    positions
+    ds.cubes[:s1_position] = YAXArray((ds.time_sentinel_1_rtc,),positions)
+    nothing
+end
+
+_to_db(x) = 10.0 * log10(x)
+
+function arceme_radar_db(ds)
+    ds.cubes[:vv_db] = _to_db.(ds.vv)
+    ds.cubes[:vh_db] = _to_db.(ds.vh)
+end
\ No newline at end of file
diff --git a/src/spatialdebias.jl b/src/spatialdebias.jl
index a480847..365c584 100644
--- a/src/spatialdebias.jl
+++ b/src/spatialdebias.jl
@@ -32,7 +32,7 @@ Computes a cloud-biased corrected footprint aggregated by land cover class for t
 provided inputcube, cloud mask and land cover cube. 
 """
 function arceme_bias_corrected_fp(band, dataset, timeaxis=:time_sentinel_2_l2a)
-    lccube = dataset.ESA_LC[time=1]
+    lccube = lckeymap.(dataset.ESA_LC[time=1])
     cloudcube = dataset.cloud_mask
     sclcube = dataset.SCL
     inputcube = dataset[band]
@@ -43,38 +43,34 @@ function arceme_bias_corrected_fp(band, dataset, timeaxis=:time_sentinel_2_l2a)
     clearsky_expected = xmap(pars ⊘ :param, fitmat ⊘ :param,inplace=false) do p,t
         p[1]*t[1]+p[2]*t[2]+p[3]*t[3]
     end
-    win_clearsky = YAXArrays.windows(clearsky_expected,lccube,expected_groups=0:100)
+    win_clearsky = YAXArrays.windows(clearsky_expected,lccube,expected_groups=1:11)
     fp_clearsky_expected = mean.(win_clearsky)[:,1,1,1,:].data
 
     clouded_expected = xmap(pars ⊘ :param, fitmat ⊘ :param, cloudcube, sclcube, inplace=false) do p, t, cl, scl
-        if (cl > 0 || scl in (1, 3, 7, 8, 9, 10, 11))
+        if _is_cloud(cl,scl)
             return NaN
         else
             p[1]*t[1]+p[2]*t[2]+p[3]*t[3]
         end
     end
-    win_clouded = YAXArrays.windows(clouded_expected,lccube,expected_groups=0:100)
+    win_clouded = YAXArrays.windows(clouded_expected,lccube,expected_groups=1:11)
     fp_clouded_expected = mean.(win_clouded)[:,1,1,1,:].data
 
     inputcube_filtered = xmap(inputcube, cloudcube,sclcube,inplace=false,output=XOutput(outtype=Float32)) do x,cl,scl
-        (cl > 0 || scl in (1, 3, 7, 8, 9, 10, 11)) ? NaN : x
+        _is_cloud(cl,scl) ? NaN : x
     end
-    win_data = YAXArrays.windows(inputcube_filtered,lccube,expected_groups=0:100)
+    win_data = YAXArrays.windows(inputcube_filtered,lccube,expected_groups=1:11)
     fp = mean.(win_data)[:,1,1,:].data
 
-    classkeys = collect(keys(arceme_classes))[2:end]
-    abundance = counter(lccube)
-    sort!(classkeys,by=i->(abundance[i],i),rev=true)
-    newdata = fp[classkeys,:] .+ fp_clearsky_expected[classkeys,:] .- fp_clouded_expected[classkeys,:]
-    classax = DD.Dim{:lc}([arceme_classes[i] for i in classkeys])
+    newdata = fp[:,:] .+ fp_clearsky_expected[:,:] .- fp_clouded_expected[:,:]
+    classax = DD.Dim{:lc}(collect(values(arceme_classes))[2:end])
 
     Dataset(
         fp = YAXArray((classax,timdim),newdata),
-        fp_uncorrected=YAXArray((classax, timdim), fp[classkeys, :]),
-        fp_clearsky_expected=YAXArray((classax, timdim), fp_clearsky_expected[classkeys, :]),
-        fp_clouded_expected=YAXArray((classax, timdim), fp_clouded_expected[classkeys, :]),
+        fp_uncorrected=YAXArray((classax, timdim), fp[:, :]),
+        fp_clearsky_expected=YAXArray((classax, timdim), fp_clearsky_expected[:, :]),
+        fp_clouded_expected=YAXArray((classax, timdim), fp_clouded_expected[:, :]),
         params=pars,
-        lc_fractions=YAXArray((classax,), [abundance[i] for i in classkeys]),
         smooth_matrix=fitmat,
     )
     
@@ -87,11 +83,11 @@ Computes a cloud-biased corrected footprint aggregated by land cover class for t
 provided inputcube, cloud mask and land cover cube. 
 """
 function arceme_uncorrected_fp(band, dataset;timeaxis=:time_sentinel_1_rtc)
-    lccube = dataset.ESA_LC[time=1]
+    lccube = lckeymap.(dataset.ESA_LC[time=1])
     inputcube = dataset[band]
     timdim = DD.dims(inputcube,timeaxis)
     
-    win = YAXArrays.windows(inputcube,lccube,expected_groups=0:100)
+    win = YAXArrays.windows(inputcube,lccube,expected_groups=1:11)
     fp = mean.(win)[:,1,1,:].data
 
     classkeys = collect(keys(arceme_classes))[2:end]
diff --git a/src/spectral_helpers.jl b/src/spectral_helpers.jl
index fe072f0..ad0c677 100644
--- a/src/spectral_helpers.jl
+++ b/src/spectral_helpers.jl
@@ -1,21 +1,37 @@
 import SpectralIndices as SI
-normalize_s2name(n) = match(r"B\d",n) === nothing ? Symbol(n) : Symbol("B0$(last(n))")
+import DiskArrayEngine as DAE
+import OnlineStats
+import CSV
+
+function load_sigmadata()
+    CSV.read(joinpath(@__DIR__, "..", "data", "sigmalist.csv"), NamedTuple)
+end
+
+const sigmadata = load_sigmadata()
+
+
+normalize_name(n) = match(r"B\d", n) !== nothing ? Symbol("B0$(last(n))") :
+                    in(n, ("VH", "VV")) ? Symbol(lowercase(n)) : Symbol(n)
 struct NTWrapper{F,names,C,I} <: Function
     f::F
     names::Val{names}
     consts::C
     indices::I
 end
-function (ntw::NTWrapper{F,names})(cl,scl,args...) where {F,names} 
-    ntw.f(cl,scl,NamedTuple{names}(args),ntw.consts,ntw.indices) 
+function (ntw::NTWrapper{F,names})(args...) where {F,names}
+    ntw.f(NamedTuple{names}(args), ntw.consts, ntw.indices)
 end
 NTWrapper(f,names::Tuple,consts::NamedTuple,indices) = NTWrapper(f,Val(names),consts,indices)
-function inner_compute_indices(cl,scl,bands,consts,indices_tuple)
-    (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return map(_ -> NaN, indices_tuple)
+function inner_compute_indices_s2(bands, consts, indices_tuple)
+    (bands.cl > 0 || (bands.scl in (1, 3, 7, 8, 9, 10, 11))) && return map(_ -> NaN, indices_tuple)
     bandparams = map(ARCEMEAnalysis.boa,bands)
     allparams = (; bandparams..., consts...)
     listofindices(indices_tuple, allparams)
 end
+function inner_compute_indices_s1(bands, consts, indices_tuple)
+    allparams = (; bands..., consts...)
+    listofindices(indices_tuple, allparams)
+end
 
 function _compute_indices(ds,indices,platform)
     indices_tuple = ((SI.indices[k] for k in indices)...,)
@@ -26,13 +42,79 @@ function _compute_indices(ds,indices,platform)
     needed_constants = intersect(allvars,keys(SpectralIndices.constants))
     undefined_constants = setdiff(allvars,needed_bands, needed_constants)
     default_values = (;(Symbol(c)=>SI.constants[c].default for c in needed_constants)...)
-    bands_to_load = (sort([normalize_s2name(SI.bands[b].platforms[platform].band) for b in needed_bands])...,)
+    bands_to_load = (sort([normalize_name(SI.bands[b].platforms[platform].band) for b in needed_bands])...,)
+    bandargs = map(b -> ds.cubes[b], bands_to_load)
+    output = map(_ -> XOutput(; outtype=Float32), indices_tuple)
+    if startswith(platform, "sentinel2")
+        wrapped_function = NTWrapper(inner_compute_indices_s2, (:cl, :scl, Symbol.(needed_bands)...), default_values, indices_tuple)
+        NamedTuple{(Symbol.(indices)...,)}(xmap(wrapped_function, ds.cloud_mask, ds.SCL, bandargs...; output, inplace=false))
+    elseif startswith(platform, "sentinel1")
+        wrapped_function = NTWrapper(inner_compute_indices_s1, (Symbol.(needed_bands)...,), default_values, indices_tuple)
+        res = xmap(wrapped_function, bandargs...; output, inplace=false)
+        length(indices) == 1 && (res = (res,))
+        NamedTuple{(Symbol.(indices)...,)}(res)
+    end
+end
 
-    wrapped_function = NTWrapper(inner_compute_indices, (Symbol.(needed_bands)...,),default_values,indices_tuple)
+#Helper functions to compute the indices from named tuples for type stability
+compute_indexx(index::SpectralIndices.SpectralIndex{<:Any,B}, params::NamedTuple) where B = index.compute(Float64, params[B]...)
+function listofindices(indices, values)
+    map(indices) do index
+        compute_indexx(index, values)
+    end
+end
 
-    bandargs = map(b->ds.cubes[b],bands_to_load)
-    output = map(_ -> XOutput(; outtype=Float32), indices_tuple)
+#Inner function to compute absolute difference of N and R band
+function absnormdiff(b1, b2, cl, scl)
+    if _is_cloud(cl, scl)
+        NaN
+    else
+        abs(boa(b1) - boa(b2))
+    end
+end
 
 
-    return NamedTuple{(Symbol.(indices)...,)}(xmap(wrapped_function, ds.cloud_mask, ds.SCL, bandargs...; output, inplace=false))
-end
\ No newline at end of file
+function estimate_sigma_per_lc(ds;bands=(:B08,:B04))
+    stat = DAE.DerivedOnlineStat{OnlineStats.ExpandingHist,OnlineStats.median,OnlineStats.fit!,(200,)}
+    f = DAE.disk_onlinestat(stat, Union{Float64,Missing}, identity, ARCEMEAnalysis.lckeymap)
+
+    si = size(ds[first(bands)].data)
+    iw1 = DAE.InputArray(absnormdiff.(ds[first(bands)].data, ds[last(bands)].data, ds.cloud_mask.data, ds.SCL.data))
+    iw2 = DAE.InputArray(ds.ESA_LC[time=1].data)
+    ow = DAE.create_outwindows((si..., 12), windows=(fill.(1, si)..., [1:12]))
+    op = DAE.GMDWop((iw1, iw2), (ow,), f)
+    res = DAE.compute(DAE.results_as_diskarrays(op)[1])
+    res[1, 1, 1, :]
+end
+
+
+function _kNDVI(cl, scl, nir, red, lc, sigmadata)
+    _is_cloud(cl, scl) && return NaN
+    sigma = sigmadata[lckeymap(lc)]
+    banddiff = boa(nir) - boa(red)
+    return tanh(((banddiff) / (2 * sigma))^2) * sign(banddiff)
+end
+
+"""
+
+"""
+function arceme_kndvi(ds, sigmadata=estimate_sigma_per_lc(ds))
+    kndvi = xmap(_kNDVI, ds.cloud_mask, ds.SCL, ds.B08, ds.B04, ds.ESA_LC[time=1], output=XOutput(outtype=Float32), function_args=(sigmadata,), inplace=false)
+    ds.cubes[:kndvi] = kndvi
+end
+
+
+
+function arceme_kndvi_pair(ds1, ds2)
+    missmean(x1, x2) =
+        if ismissing(x1)
+            ismissing(x2) ? missing : x2
+        else
+            ismissing(x2) ? x1 : 0.5 * (x1 + x2)
+        end
+    sigmadata1 = estimate_sigma_per_lc(ds1)
+    sigmadata2 = estimate_sigma_per_lc(ds2)
+    sigmadata = missmean.(sigmadata1, sigmadata2)
+    arceme_kndvi(ds1, sigmadata)
+    arceme_kndvi(ds2, sigmadata)
+end

From 52485117f2f79239d9e7108d22604eac5300a306 Mon Sep 17 00:00:00 2001
From: Fabian Gans <fgans@bgc-jena.mpg.de>
Date: Mon, 23 Feb 2026 12:09:27 +0000
Subject: [PATCH 8/8] Remove precomputed sigmas

---
 src/ARCEMEAnalysis.jl   | 2 +-
 src/spectral_helpers.jl | 7 -------
 2 files changed, 1 insertion(+), 8 deletions(-)

diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index b942e9b..d485324 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -417,7 +417,7 @@ function arceme_create_indexcubes(; indices_s1=["DpRVIVV"], indices_s2=["NDVI",
             output_base = "$local_cubepath/ARCEME-DC-6-INDICES"
             name = arceme_cubename(event)
             indexcube = setchunks(ds[fields_to_save], (500, 500, 25))
-            compute_to_zarr(indexcube, joinpath(output_base, name), overwrite=true)
+            compute_to_zarr(indexcube, joinpath(output_base, name), custom_loopranges=(500, 500, 25), overwrite=true)
             cube2 = setchunks(ds[["vv_db", "vh_db", "kNDVI"]], (500, 500, 25))
             savedataset(cube2, path=joinpath(output_base, name), append=true)
             run(Cmd(`zip -0 -r ../$(name).zip .`, dir=joinpath(output_base, name)))
diff --git a/src/spectral_helpers.jl b/src/spectral_helpers.jl
index ad0c677..7c3fd28 100644
--- a/src/spectral_helpers.jl
+++ b/src/spectral_helpers.jl
@@ -3,13 +3,6 @@ import DiskArrayEngine as DAE
 import OnlineStats
 import CSV
 
-function load_sigmadata()
-    CSV.read(joinpath(@__DIR__, "..", "data", "sigmalist.csv"), NamedTuple)
-end
-
-const sigmadata = load_sigmadata()
-
-
 normalize_name(n) = match(r"B\d", n) !== nothing ? Symbol("B0$(last(n))") :
                     in(n, ("VH", "VV")) ? Symbol(lowercase(n)) : Symbol(n)
 struct NTWrapper{F,names,C,I} <: Function