diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..db7bc27
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1 @@
+style = "default"
\ No newline at end of file
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 81bf6f5..833e3ae 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -34,7 +34,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
   docs:
@@ -49,7 +49,7 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - name: Configure doc environment
         shell: julia --project=docs --color=yes {0}
         run: |
diff --git a/Project.toml b/Project.toml
index d9fbaa4..fff9fe6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,14 +8,20 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DateFormats = "44557152-fe0a-4de1-8405-416d90313ce6"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
-GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
+DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Minio = "4281f0d9-7ae0-406e-9172-b7277c1efa20"
@@ -24,7 +30,9 @@ Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+RollingFunctions = "b0e4dd01-7b14-53d8-9b45-175a3e362653"
 SpectralIndices = "df0093a1-273d-40bc-819a-796ec3476907"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StringViews = "354b36f9-a18e-4713-926e-db85100087ba"
@@ -37,14 +45,20 @@ ArchGDAL = "0.10.11"
 CSV = "0.10.15"
 Colors = "0.13.1"
 CondaPkg = "0.2.34"
+DataFrames = "1.8.2"
 DataStructures = "0.18.22, 0.19"
+DateFormats = "0.1.20"
 Dates = "1.11.0"
 Dierckx = "0.5.4"
 DimensionalData = "0.29.24"
 DiskArrayEngine = "0.3.0"
-GeoMakie = "0.7.16"
+DiskArrays = "0.4.21"
+Distributed = "1.11.0"
 GeoJSON = "0.8.4"
+GeoMakie = "0.7.16"
 HTTP = "1.10.19"
+HypothesisTests = "0.11.7"
+Interpolations = "0.15.1"
 JSON = "1.4.0"
 Makie = "0.24.9"
 Minio = "0.2.2"
@@ -53,11 +67,13 @@ Preferences = "1.5.1"
 ProgressMeter = "1.11.0"
 Proj = "1.9.0"
 PythonCall = "0.9.31"
+Random = "1.11.0"
 Reexport = "1.2.2"
+RollingFunctions = "0.8.1"
 SpectralIndices = "0.2.15"
 Statistics = "1.11.1"
 StringViews = "1.3.7"
-YAXArrays = "0.7.1"
+YAXArrays = "0.7.7"
 Zarr = "0.9.5"
 ZipArchives = "2.6.0"
 julia = "1.10"
diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index d8ab65c..b45e995 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -1,23 +1,40 @@
 module ARCEMEAnalysis
 using Zarr: S3Store, Zarr
 using Minio: MinioConfig
-using YAXArrays: open_dataset, ⊘, xmap, XOutput, YAXArray, YAXArrays, Dataset, setchunks, compute_to_zarr, savedataset
+using YAXArrays:
+    open_dataset,
+    ⊘,
+    xmap,
+    XOutput,
+    YAXArray,
+    YAXArrays,
+    Dataset,
+    setchunks,
+    compute_to_zarr,
+    savedataset
 import DimensionalData as DD
 using Colors: RGB, RGBA
 using Dates: DateTime, Year, Date
 import CSV
-using Statistics: mean
+using Statistics: mean, median
 using DataStructures: SortedDict, counter
 using ProgressMeter: @showprogress
 using SpectralIndices: compute_index, SpectralIndices
 import Proj
 import SpectralIndices as SI
 import GeoJSON
+import Distributed: nprocs, pmap
 
 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))
+    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)
@@ -25,14 +42,10 @@ end
 
 const arceme_legends = Dict(
     # Major Main Crop Type
-    "MCTY" => SortedDict(
-        1 => "No cropland",
-        2 => "Annual crop",
-        3 => "Permanent crop"
-    ),
+    "MCTY" => SortedDict(1 => "No cropland", 2 => "Annual crop", 3 => "Permanent crop"),
     # Main Crop Type
     "CTY" => SortedDict(
-        0 =>  "No cropland",
+        0 => "No cropland",
         1110 => "Wheat",
         1120 => "Barley",
         1130 => "Maize",
@@ -52,40 +65,37 @@ const arceme_legends = Dict(
         2320 => "Nuts",
         3100 => "Unclassified annual crop",
         3200 => "Unclassified permanent crop",
-        65535 => "Outside area"
+        65535 => "Outside area",
     ),
     # Secondary Crop Type
     "CPSCT" => SortedDict(
-        0 => "No annual cropland", 
+        0 => "No annual cropland",
         1 => "Short Summer",
         2 => "Long Summer",
         3 => "Short Winter",
-        4 => "Long Winter", 
-        65526 => "Fallow land", 
-        65527 => "No cropping pattern detected",  
+        4 => "Long Winter",
+        65526 => "Fallow land",
+        65527 => "No cropping pattern detected",
         65530 => "No secondary crop growing season delineated",
-        65531 => "Not enough data", 
-        65532 => "No cropping season detected",  
+        65531 => "Not enough data",
+        65532 => "No cropping season detected",
         65533 => "Growing season extends beyond timeframe",
-        65535 => "Outside area"
+        65535 => "Outside area",
     ),
     # Fallow Land
-    "CPFLP" => SortedDict(
-        0 => "No fallow land",
-        1 => "Fallow land",
-        65535 => "Outside area"
-    ),
+    "CPFLP" =>
+        SortedDict(0 => "No fallow land", 1 => "Fallow land", 65535 => "Outside area"),
     # Cropping Seasons in Year
     "CPCSY" => SortedDict(
-        0 => "No annual cropland", 
+        0 => "No annual cropland",
         1 => "One growing season",
         2 => "Two growing seasons",
-        65526 => "Fallow land", 
-        65527 => "No cropping pattern detected",  
-        65531 => "Not enough data", 
-        65532 => "No cropping season detected",  
+        65526 => "Fallow land",
+        65527 => "No cropping pattern detected",
+        65531 => "Not enough data",
+        65532 => "No cropping season detected",
         65533 => "Growing season extends beyond timeframe",
-        65535 => "Outside area"
+        65535 => "Outside area",
     ),
     # Tree cover: Leaf Type
     "DLT" => SortedDict(
@@ -97,46 +107,67 @@ const arceme_legends = Dict(
     # Grassland
     "GRA" => SortedDict(
         0 => "all non-grassland areas",
-        1 => "grassland", 
+        1 => "grassland",
         255 => "outside area",
     ),
     "ESA_LC" => SortedDict(
-    0 => "No data",
-    10 => "Tree cover",
-    20 => "Shrubland",
-    30 => "Grassland",
-    40 => "Cropland",
-    50 => "Built-up",
-    60 => "Bare /sparse vegetation",
-    70 => "Snow and Ice",
-    80 => "Permanent water bodies",
-    90 => "Herbaceous wetland",
-    95 => "Mangroves",
-    100 => "Moss and lichen",
-    )
+        0 => "No data",
+        10 => "Tree cover",
+        20 => "Shrubland",
+        30 => "Grassland",
+        40 => "Cropland",
+        50 => "Built-up",
+        60 => "Bare /sparse vegetation",
+        70 => "Snow and Ice",
+        80 => "Permanent water bodies",
+        90 => "Herbaceous wetland",
+        95 => "Mangroves",
+        100 => "Moss and lichen",
+    ),
 )
-
-lckeymap(k::Any) = ifelse(k > 90, (k + 10), k) ÷ 10 + 1
-function lckeymap(ds::Dataset;strata="ESA_LC")
-    if strata == "ESA_LC" && return lckeymap.(ds.ESA_LC[time=1]) end
-    if strata == "CTY" && return HRL.ctykeymap.(ds.CTY) end
-    if strata == "MCTY" && return HRL.mctykeymap.(ds.CTY) end
+lckeymap(k) = ifelse(k > 90, (k + 10), k) ÷ 10 + 1
+function lckeymap(ds::Dataset; strata = "ESA_LC")
+    if strata == "ESA_LC" && return lckeymap.(ds.ESA_LC[time=1])
+    end
+    if strata == "CTY" && return HRL.ctykeymap.(ds.CTY)
+    end
+    if strata == "MCTY" && return HRL.mctykeymap.(ds.CTY)
+    end
 end
 
 
-export arceme_cubename, arceme_open, arceme_starttime, arceme_endtime, arceme_eventdate,
-    arceme_coordinates, arceme_ndvi, arceme_rgb, arceme_eventlist, arceme_eventpairs,
-    arceme_landcover, arceme_optical_band_fingerprints, arceme_radar_fingerprints,
-    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral, arceme_kndvi, arceme_radar_db,
-    arceme_create_indexcubes, arceme_index_fingerprints, arceme_open_fingerprint, arceme_fractions,
-    arceme_legends
+export arceme_cubename,
+    arceme_open,
+    arceme_starttime,
+    arceme_endtime,
+    arceme_eventdate,
+    arceme_coordinates,
+    arceme_ndvi,
+    arceme_rgb,
+    arceme_eventlist,
+    arceme_eventpairs,
+    arceme_landcover,
+    arceme_optical_band_fingerprints,
+    arceme_radar_fingerprints,
+    time_aggregate_fingerprint,
+    arceme_validpairs,
+    arceme_spectral,
+    arceme_kndvi,
+    arceme_radar_db,
+    arceme_create_indexcubes,
+    arceme_index_fingerprints,
+    arceme_open_fingerprint,
+    arceme_fractions,
+    arceme_legends,
+    arceme_cloudmask_halo,
+    arceme_precompute_fingerprints
 
 """
     _arceme_cubenames(;batch="6")
 
 List all available data cube names in the specified batch stored in the ARCEME S3 bucket.
 """
-function _arceme_cubenames(; batch="ARCEME-DC-6")
+function _arceme_cubenames(; batch = "ARCEME-DC-DHP-GLOBAL")
     cubenames = if local_cubepath === nothing
         store = S3Store("$(batch)/", MinioConfig(httpstore))
 
@@ -180,7 +211,9 @@ end
 
 Get a list of ARCEME events from the local CSV or geojson file.
 """
-function arceme_eventlist(;csv=joinpath(@__DIR__, "..","data","dhp_global_subselection.csv"))
+function arceme_eventlist(;
+    csv = joinpath(@__DIR__, "..", "data", "dhp_global_subselection.csv"),
+)
     map(CSV.File(csv)) do row
         Event(
             row.uid,
@@ -200,7 +233,10 @@ function arceme_eventlist(;csv=joinpath(@__DIR__, "..","data","dhp_global_subsel
     end
 end
 
-function arceme_eventlist(source::String; geojson=joinpath(@__DIR__, "..","data/selection_eu_wocat_dhp_qdoy_hrl.geojson")) 
+function arceme_eventlist(
+    source::String;
+    geojson = joinpath(@__DIR__, "..", "data/selection_eu_wocat_dhp_qdoy_hrl.geojson"),
+)
     map(GeoJSON.read(geojson)) do row
         Event(
             row.uid,
@@ -209,13 +245,12 @@ function arceme_eventlist(source::String; geojson=joinpath(@__DIR__, "..","data/
             row.geometry[2],
             DateTime(row.startdate),
             Symbol(source),
-
         )
     end
 end
 
-arceme_cubename(event::Event) = 
-"DC__$(event.uid)__$(Date(arceme_starttime(event)))__$(Date(arceme_endtime(event))).zarr"
+arceme_cubename(event::Event) =
+    "DC__$(event.uid)__$(Date(arceme_starttime(event)))__$(Date(arceme_endtime(event))).zarr"
 
 """
     arceme_eventpairs()
@@ -223,33 +258,65 @@ arceme_cubename(event::Event) =
 Get pairs of ARCEME events from the event list.
 """
 arceme_eventpairs() =
-    Iterators.partition(sort(arceme_eventlist(),by=i->(i.event_label,i.source)),2) |> collect
+    Iterators.partition(sort(arceme_eventlist(), by = i -> (i.event_label, i.source)), 2) |>
+    collect
 
 """
-    arceme_validpairs(;batch="ARCEME-DC-6")
+    arceme_validpairs(;batch="ARCEME-DC-DHP-GLOBAL")
 
 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")
+function arceme_validpairs(; batch = "ARCEME-DC-DHP-GLOBAL")
     allpairs = arceme_eventpairs()
     validpairs = if local_cubepath === nothing
-        map(x -> all([Zarr.is_zgroup(Zarr.HTTPStore("$httpstore/$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)
+        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; batch="ARCEME-DC-6", indices=true, trylocal=true, fingerprint=true, hrl=true)
-    arceme_open(arceme_cubename(event); batch=batch, indices=indices, trylocal=trylocal, fingerprint=fingerprint, hrl=true)
+function arceme_open(
+    event::Event;
+    batch = "ARCEME-DC-DHP-GLOBAL",
+    indices = true,
+    trylocal = true,
+    fingerprint = true,
+    hrl = true,
+)
+    arceme_open(
+        arceme_cubename(event);
+        batch = batch,
+        indices = indices,
+        trylocal = trylocal,
+        fingerprint = fingerprint,
+        hrl = true,
+    )
 end
 
-function arceme_open_fingerprint(event::Event; batch="ARCEME-DC-6")
+function arceme_open_fingerprint(event::Event; batch = "ARCEME-DC-DHP-GLOBAL")
     arceme_open_fingerprint(arceme_cubename(event); batch)
 end
 
-arceme_open_fingerprint(cubename; batch="ARCEME-DC-6") = open_dataset(joinpath(local_cubepath, "$batch-fingerprints.zip"), path=cubename)
+arceme_open_fingerprint(cubename; batch = "ARCEME-DC-DHP-GLOBAL") = open_dataset(
+    joinpath(local_cubepath, "$batch-fingerprints.zip"),
+    path = cubename,
+    force_datetime = true,
+)
 
 """
      arceme_landcover(ds)
@@ -260,15 +327,21 @@ function arceme_landcover(ds)
     cdr = counter(ds.ESA_LC)
     map(collect(arceme_legends["ESA_LC"])) do (k, v)
         count = get(cdr, k, 0)
-        (key=k, class=v, count=count, fraction=count / 1000000)
+        (
+            key = k,
+            class = v,
+            count = count,
+            fraction = count / (length(ds.x) * length(ds.y)),
+        )
     end
 end
-arceme_landcover(ev::Event; batch="ARCEME-DC-6") = arceme_landcover(arceme_open(ev, batch=batch))
+arceme_landcover(ev::Event; batch = "ARCEME-DC-DHP-GLOBAL") =
+    arceme_landcover(arceme_open(ev, batch = batch))
 
 
 """ 
-    arceme_open(cubename::String; batch="ARCEME-DC-6", indices=true, hrl=true)
-    arceme_open(event::Event; batch="ARCEME-DC-6", indices=true, hrl=true)
+    arceme_open(cubename::String; batch="ARCEME-DC-DHP-GLOBAL", indices=true, hrl=true)
+    arceme_open(event::Event; batch="ARCEME-DC-DHP-GLOBAL", indices=true, hrl=true)
 
 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). 
@@ -276,13 +349,59 @@ If `indices=true` (`hrl=true`), also opens indices (hrl layers) if they exist on
 
 See also `arceme_create_indexcubes` and `HRL.hrl_warp`.
 """
-function arceme_open(cubename; batch="ARCEME-DC-6", indices=true, trylocal=true, fingerprint=true, hrl=true)
+function arceme_open(
+    cubename;
+    batch = "ARCEME-DC-DHP-GLOBAL",
+    indices = true,
+    trylocal = true,
+    fingerprint = true,
+    hrl = true,
+)
     if local_cubepath === nothing || !trylocal
-        open_dataset("$httpstore/$batch/$cubename", force_datetime=true)
+        open_dataset("$httpstore/$batch/$cubename", force_datetime = true)
     else
-        main_ds = open_dataset(joinpath(local_cubepath, batch, string(cubename, ".zip")))
-        if indices && isfile(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
-            index_ds = open_dataset(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
+        main_ds = open_dataset(
+            joinpath(local_cubepath, batch, string(cubename, ".zip")),
+            force_datetime = true,
+        )
+        if batch == "ARCEME-DC-DHP-GLOBAL"
+            # if :ESA_LC ∉ keys(main_ds.cubes)
+            #     @show cubename
+            #     oldds = open_dataset(joinpath(local_cubepath, "ARCEME-DC-6", string(cubename, ".zip")), force_datetime=true)
+            #     main_ds.cubes[:ESA_LC] = oldds.ESA_LC
+            # end
+            if :B01 ∉ keys(main_ds.cubes)
+                oldds = open_dataset(
+                    joinpath(local_cubepath, "ARCEME-DC-6", string(cubename, ".zip")),
+                    force_datetime = true,
+                )
+                main_ds.axes[:time_sentinel_2_l2a] = oldds.axes[:time_sentinel_2_l2a]
+                for v in [
+                    :B01,
+                    :B02,
+                    :B03,
+                    :B04,
+                    :B05,
+                    :B06,
+                    :B07,
+                    :B08,
+                    :B09,
+                    :B11,
+                    :B12,
+                    :B8A,
+                    :SCL,
+                    :cloud_mask,
+                ]
+                    main_ds.cubes[v] = oldds.cubes[v]
+                end
+            end
+        end
+        if indices &&
+           isfile(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
+            index_ds = open_dataset(
+                joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")),
+                force_datetime = true,
+            )
             for (k, v) in (index_ds.cubes)
                 main_ds.cubes[k] = v
             end
@@ -290,13 +409,19 @@ function arceme_open(cubename; batch="ARCEME-DC-6", indices=true, trylocal=true,
         if fingerprint
             try
                 ds_fp = arceme_open_fingerprint(cubename; batch)
-                main_ds = Dataset(; main_ds.cubes..., ds_fp.cubes...)
+                main_ds = Dataset(;
+                    main_ds.cubes...,
+                    ds_fp.cubes...,
+                    properties = main_ds.properties,
+                )
             catch
                 main_ds
             end
         end
         if hrl && isfile(joinpath(local_cubepath, "$batch-HRL", string(cubename, ".zip")))
-            hrl_ds = open_dataset(joinpath(local_cubepath, "$batch-HRL", string(cubename, ".zip")))
+            hrl_ds = open_dataset(
+                joinpath(local_cubepath, "$batch-HRL", string(cubename, ".zip")),
+            )
             for (k, v) in (hrl_ds.cubes)
                 main_ds.cubes[k] = v
             end
@@ -339,7 +464,7 @@ Get the longitude and latitude coordinates of the center of the ARCEME data cube
 function arceme_coordinates(ds)
     trans = Proj.Transformation("EPSG:$(ds.properties["epsg"])", "OGC:84")
     x, y = ds.properties["central_x"], ds.properties["central_y"]
-    lon, lat = round.(trans((x, y)), digits=4)
+    lon, lat = round.(trans((x, y)), digits = 4)
     return lon, lat
 end
 
@@ -365,7 +490,7 @@ end
 
 Compute the listed indices using SpectralIndices.jl. Not working.
 """
-function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2")
+function arceme_spectral(ds, indices::Vector{String}; platform = "sentinel2")
     tmp = if platform == "sentinel2" || platform == "sentinel2a" || platform == "sentinel2b"
         pl = platform == "sentinel2" ? "sentinel2a" : platform
         _compute_indices(ds, indices, pl)
@@ -381,7 +506,7 @@ function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2")
 end
 
 _is_cloud(cl, scl) = (cl > 0 || (scl in (0, 1, 3, 7, 8, 9, 10, 11)))
-
+_is_clear_or_shadow(cl, scl) = !(cl in (1, 2) || (scl in (0, 1, 3, 7, 8, 9, 10, 11)))
 
 
 """
@@ -391,20 +516,20 @@ Compute the radiometric offsets for Sentinel 2 bands to get values at bottom of
 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
+boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000) =
+    (band + BOA_ADD_OFFSET) / QUANTIFICATION_VALUE
 
 """
     arceme_rgb(ds)
 
 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)
-        # 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(clamp(boa(r), 0, 1), clamp(boa(g), 0, 1), clamp(boa(b), 0, 1))
-    end
+arceme_rgb(ds) = broadcast(ds.B02, ds.B03, ds.B04) do b, g, r
+    # 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(clamp(boa(r), 0, 1), clamp(boa(g), 0, 1), clamp(boa(b), 0, 1))
+end
 
 """
     arceme_optical_band_fingerprints(ds_d, ds_dhp)
@@ -423,18 +548,37 @@ function arceme_optical_band_fingerprints(ds_d, ds_dhp)
         arceme_bias_corrected_fp(band, ds_d)
     end
     fingerprint_sparse_d = YAXArrays.concatenatecubes(map(i -> i.fp, fp_d), banddim)
-    fingerprint_uncor_sparse_d = YAXArrays.concatenatecubes(map(i -> i.fp_uncorrected, fp_d), banddim)
+    fingerprint_uncor_sparse_d =
+        YAXArrays.concatenatecubes(map(i -> i.fp_uncorrected, fp_d), banddim)
 
     fp_dhp = @showprogress desc = "DHP fingerprint......" map(optical_bands) do band
         arceme_bias_corrected_fp(band, ds_dhp)
     end
     fingerprint_sparse_dhp = YAXArrays.concatenatecubes(map(i -> i.fp, fp_dhp), banddim)
-    fingerprint_uncor_sparse_dhp = YAXArrays.concatenatecubes(map(i -> i.fp_uncorrected, fp_dhp), banddim)
-
-    fingerprints_d = time_aggregate_fingerprint(fingerprint_sparse_d, eventdate_d, banddim, :time_sentinel_2_l2a)
-    fingerprints_dhp = time_aggregate_fingerprint(fingerprint_sparse_dhp, eventdate_dhp, banddim, :time_sentinel_2_l2a)
+    fingerprint_uncor_sparse_dhp =
+        YAXArrays.concatenatecubes(map(i -> i.fp_uncorrected, fp_dhp), banddim)
+
+    fingerprints_d = time_aggregate_fingerprint(
+        fingerprint_sparse_d,
+        eventdate_d,
+        banddim,
+        :time_sentinel_2_l2a,
+    )
+    fingerprints_dhp = time_aggregate_fingerprint(
+        fingerprint_sparse_dhp,
+        eventdate_dhp,
+        banddim,
+        :time_sentinel_2_l2a,
+    )
 
-    Dataset(; fingerprints_d, fingerprints_dhp, fingerprint_sparse_d, fingerprint_sparse_dhp, fingerprint_uncor_sparse_d, fingerprint_uncor_sparse_dhp)
+    Dataset(;
+        fingerprints_d,
+        fingerprints_dhp,
+        fingerprint_sparse_d,
+        fingerprint_sparse_dhp,
+        fingerprint_uncor_sparse_d,
+        fingerprint_uncor_sparse_dhp,
+    )
 end
 
 """
@@ -461,10 +605,25 @@ function arceme_radar_fingerprints(ds_d, ds_dhp)
     end
     fingerprint_sparse_dhp = YAXArrays.concatenatecubes(map(i -> i.fp, fp_dhp), banddim)
 
-    fingerprints_d = time_aggregate_fingerprint(fingerprint_sparse_d, eventdate_d, banddim, :time_sentinel_1_rtc)
-    fingerprints_dhp = time_aggregate_fingerprint(fingerprint_sparse_dhp, eventdate_dhp, banddim, :time_sentinel_1_rtc)
+    fingerprints_d = time_aggregate_fingerprint(
+        fingerprint_sparse_d,
+        eventdate_d,
+        banddim,
+        :time_sentinel_1_rtc,
+    )
+    fingerprints_dhp = time_aggregate_fingerprint(
+        fingerprint_sparse_dhp,
+        eventdate_dhp,
+        banddim,
+        :time_sentinel_1_rtc,
+    )
 
-    Dataset(; fingerprints_d, fingerprints_dhp, fingerprint_sparse_d, fingerprint_sparse_dhp)
+    Dataset(;
+        fingerprints_d,
+        fingerprints_dhp,
+        fingerprint_sparse_d,
+        fingerprint_sparse_dhp,
+    )
 end
 
 
@@ -476,7 +635,11 @@ Merges the optical and radar fingerprints datasets into a single dataset.
 function arceme_merge_fingerprints(optical_fp, radar_fp)
     newbands = DD.Dim{:band}([optical_fp.band.val; radar_fp.band.val])
     newarrays = map((:fingerprints_d, :fingerprints_dhp)) do cubename
-        mergeddata = cat(optical_fp[cubename].data ./ typemax(UInt16), radar_fp[cubename].data, dims=3)
+        mergeddata = cat(
+            optical_fp[cubename].data ./ typemax(UInt16),
+            radar_fp[cubename].data,
+            dims = 3,
+        )
         cubename => YAXArray((radar_fp.lc, radar_fp.Ti, newbands), mergeddata)
     end
     Dataset(; newarrays...)
@@ -499,18 +662,29 @@ function time_aggregate_fingerprint(allbands, eventdate, banddim, timeaxis)
     step_per_year = 6
     ts = Array(DD.dims(allbands, timeaxis))
     data = allbands.data
-    timestep_to_group(t, eventdate, step_per_year) = (clamp(ceil(Int, (t - eventdate).value / (365.25 * 24 * 60 * 60 * 1000) * step_per_year), -step_per_year + 1, step_per_year))
+    timestep_to_group(t, eventdate, step_per_year) = (clamp(
+        ceil(Int, (t - eventdate).value / (365.25 * 24 * 60 * 60 * 1000) * step_per_year),
+        -step_per_year + 1,
+        step_per_year,
+    ))
     groupinds = timestep_to_group.(ts, eventdate, step_per_year)
     fingerprint_timesteps = (-step_per_year+1):(step_per_year)
     res = map(fingerprint_timesteps) do t
         ii = findall(==(t), groupinds)
         isnothing(ii) && return NaN
-        mapslices(data[:, ii, :], dims=2) do sts
+        mapslices(data[:, ii, :], dims = 2) do sts
             all(ismissing, sts) && return NaN
             mean(skipmissing(sts))
         end
     end
-    YAXArray((allbands.lc, DD.Ti((fingerprint_timesteps .- 0.5) ./ step_per_year .* 12), banddim), cat(res..., dims=2))
+    YAXArray(
+        (
+            allbands.lc,
+            DD.Ti((fingerprint_timesteps .- 0.5) ./ step_per_year .* 12),
+            banddim,
+        ),
+        cat(res..., dims = 2),
+    )
 end
 
 """
@@ -520,62 +694,107 @@ For every valid event pair creates and stores a data cubes of a list of precompu
 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"], batch="ARCEME-DC-6", subset=:)
-
-    for ev in arceme_validpairs(batch=batch)[subset]
-        ds_pair = arceme_open.(ev, indices=false, batch=batch)
-
+function arceme_create_indexcubes(;
+    indices_s1 = ["DpRVIVV"],
+    indices_s2 = ["NDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"],
+    batch = "ARCEME-DC-DHP-GLOBAL",
+    subset = :,
+)
+    @showprogress pmap(arceme_validpairs(batch = batch)[subset]) do ev
+        ds_pair = arceme_open.(ev, indices = false, batch = batch)
         foreach(ds_pair) do ds
-            arceme_spectral(ds, indices_s1, platform="sentinel1")
-            arceme_spectral(ds, indices_s2, platform="sentinel2")
-            arceme_radar_db(ds)
+            do_s1 = haskey(ds.cubes, :vh)
+            do_s1 && arceme_spectral(ds, indices_s1, platform = "sentinel1")
+            arceme_spectral(ds, indices_s2, platform = "sentinel2")
+            do_s1 && arceme_radar_db(ds)
+            arceme_cloudmask_halo(ds)
             arceme_fractions(ds)
         end
         ARCEMEAnalysis.arceme_kndvi_pair(ds_pair...)
 
-        fields_to_save = [indices_s1; indices_s2]
+
 
         output_base = "$local_cubepath/$(batch)-INDICES"
 
         foreach(ds_pair, ev) do ds, event
+            do_s1 = haskey(ds.cubes, :vh)
+            fields_to_save = do_s1 ? [indices_s1; indices_s2] : indices_s2
             writeindexcubes(ds, event, output_base, fields_to_save)
         end
     end
+    nothing
 end
 
 function writeindexcubes(ds, event, output_base, fields_to_save)
     name = arceme_cubename(event)
     indexcube = setchunks(ds[fields_to_save], (500, 500, 25))
-    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)
-    cube3 = ds[["cloud_fraction", "lc_fraction"]]
-    savedataset(cube3, path=joinpath(output_base, name), append=true)
-    isfile(joinpath(output_base, string(name, ".zip"))) && rm(joinpath(output_base, string(name, ".zip")))
-    run(Cmd(`7z a -tzip -mx=0 ../$(name).zip .`, dir=joinpath(output_base, name)))
-    rm(joinpath(output_base, name), recursive=true)
+    compute_to_zarr(
+        indexcube,
+        joinpath(output_base, name),
+        custom_loopranges = (500, 500, 25),
+        overwrite = true,
+        showprogress = false,
+        use_dagger = false,
+    )
+    if haskey(ds.cubes, :vv_db)
+        cube2 = setchunks(ds[["vv_db", "vh_db", "kNDVI"]], (500, 500, 25))
+        savedataset(cube2, path = joinpath(output_base, name), append = true)
+    end
+    cube3 = ds[["cloud_fraction", "cloudmask_halo"]]
+    savedataset(cube3, path = joinpath(output_base, name), append = true)
+    isfile(joinpath(output_base, string(name, ".zip"))) &&
+        rm(joinpath(output_base, string(name, ".zip")))
+    zip_dir(joinpath(output_base, "$name.zip"), joinpath(output_base, name))
+    rm(joinpath(output_base, name), recursive = true)
+end
+
+# Utility function to zip a directory
+function apply_dirfiles(f, basepath, dir = "")
+    flist = readdir(joinpath(basepath, dir))
+    for file in flist
+        if isdir(joinpath(basepath, dir, file))
+            apply_dirfiles(f, basepath, joinpath(dir, file))
+        else
+            f(basepath, joinpath(dir, file))
+        end
+    end
+end
+function zip_dir(output, d)
+    ZipWriter(output) do w
+        apply_dirfiles(d) do p, file
+            zip_newfile(w, file)
+            write(w, read(joinpath(p, file)))
+        end
+    end
 end
 
 """
-`arceme_create_indexcubes(event_list; indices_s1=["DpRVIVV", "vv_db", "vh_db"], indices_s2=["NDVI", "kNDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"], batch="ARCEME-DC-6")`
+`arceme_create_indexcubes(event_list; indices_s1=["DpRVIVV", "vv_db", "vh_db"], indices_s2=["NDVI", "kNDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"], batch="ARCEME-DC-DHP-GLOBAL")`
 
 Computes indices for single events (not pairs) and stores them at `joinpath(local_cubepath,\"\$(batch)-INDICES\")`. 
 """
-function arceme_create_indexcubes(event_list; indices_s1=["DpRVIVV"], indices_s2=["NDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"], batch="ARCEME-DC-6")
-    if !isdir(joinpath(local_cubepath,"$(batch)-INDICES"))
-        mkdir(joinpath(local_cubepath,"$(batch)-INDICES"))
-    end        
-    output_base = joinpath(local_cubepath,"$(batch)-INDICES")
+function arceme_create_indexcubes(
+    event_list;
+    indices_s1 = ["DpRVIVV"],
+    indices_s2 = ["NDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"],
+    batch = "ARCEME-DC-DHP-GLOBAL",
+)
+    if !isdir(joinpath(local_cubepath, "$(batch)-INDICES"))
+        mkdir(joinpath(local_cubepath, "$(batch)-INDICES"))
+    end
+    output_base = joinpath(local_cubepath, "$(batch)-INDICES")
 
     for event in event_list
         name = arceme_cubename(event)
-        isfile(joinpath(output_base, string(name, ".zip"))) && rm(joinpath(output_base, string(name, ".zip")))
+        isfile(joinpath(output_base, string(name, ".zip"))) &&
+            rm(joinpath(output_base, string(name, ".zip")))
 
-        ds = arceme_open(event, batch = batch, indices=false, hrl=false)
+        ds = arceme_open(event, batch = batch, indices = false, hrl = false)
 
-        arceme_spectral(ds, indices_s1, platform="sentinel1")
-        arceme_spectral(ds, indices_s2, platform="sentinel2")
+        arceme_spectral(ds, indices_s1, platform = "sentinel1")
+        arceme_spectral(ds, indices_s2, platform = "sentinel2")
         arceme_radar_db(ds)
+        arceme_cloudmask_halo(ds)
         arceme_fractions(ds)
         arceme_kndvi(ds)
 
@@ -586,50 +805,80 @@ function arceme_create_indexcubes(event_list; indices_s1=["DpRVIVV"], indices_s2
     end
 end
 """
-`arceme_index_fingerprints(ds; indices_s1=["DpRVIVV", "vv_db", "vh_db"], indices_s2=["NDVI", "kNDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"], strata="ESA_LC")`
+`arceme_index_fingerprints(ds; indices_s1=["DpRVIVV", "vv_db", "vh_db"], indices_s2=["NDVI", "kNDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"], strata="ESA_LC", fpstat=mean)`
 
-Computes fingerprints for indices, stratified by `strata`.
+Computes fingerprints for indices, stratified by `strata` using statistic `fpstat`.
 """
-function arceme_index_fingerprints(ds; indices_s1=["DpRVIVV", "vv_db", "vh_db"], indices_s2=["NDVI", "kNDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"], strata="ESA_LC")
+function arceme_index_fingerprints(
+    ds;
+    indices_s1 = ["DpRVIVV", "vv_db", "vh_db"],
+    indices_s2 = ["NDVI", "kNDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"],
+    strata = "ESA_LC",
+    fpstat = mean,
+)
 
     #banddim = DD.Dim{:band}(string.(optical_bands))
     indices_s2 = filter(i -> in(Symbol(i), keys(ds.cubes)), indices_s2) |> collect
     indices_s1 = filter(i -> in(Symbol(i), keys(ds.cubes)), indices_s1) |> collect
-    if isempty(indices_s1) || isempty(indices_s2)
-        return nothing
-    end
-    indexdim_s2 = DD.Dim{:band_s2}(indices_s2)
-    indexdim_s1 = DD.Dim{:band_s1}(indices_s1)
-    #eventdate = arceme_eventdate(ds)
-    fp_s2 = @showprogress desc = "S2 fingerprint.." map(indices_s2) do band
-        arceme_bias_corrected_fp(band, ds, strata=strata)
+    s1cubes, s2cubes = (), ()
+    if !isempty(indices_s1)
+        indexdim_s1 = DD.Dim{:band_s1}(indices_s1)
+        fp_s1 = map(indices_s1) do band
+            arceme_uncorrected_fp(band, ds, timeaxis = :time_sentinel_1_rtc, strata = strata)
+        end
+        fingerprint_sparse_s1 = setchunks(
+            YAXArrays.concatenatecubes(map(i -> i.fp, fp_s1), indexdim_s1),
+            (; band_s1 = length(indices_s1)),
+        )
+        s1cubes = (
+            :s1_indices => fingerprint_sparse_s1,
+            :class_fractions => fp_s1[1].class_fractions,
+        )
     end
-    fingerprint_sparse_s2 = setchunks(YAXArrays.concatenatecubes(map(i -> i.fp, fp_s2), indexdim_s2), (; band_s2=length(indices_s2)))
-    fingerprint_uncor_sparse_s2 = setchunks(YAXArrays.concatenatecubes(map(i -> i.fp_uncorrected, fp_s2), indexdim_s2), (; band_s2=length(indices_s2)))
-    #fingerprints_s2 = time_aggregate_fingerprint(fingerprint_sparse_s2, eventdate, indexdim_s2, :time_sentinel_2_l2a)
+    if !isempty(indices_s2)
+        indexdim_s2 = DD.Dim{:band_s2}(indices_s2)
 
-    fp_s1 = @showprogress desc = "S1 fingerprint.." map(indices_s1) do band
-        arceme_uncorrected_fp(band, ds, timeaxis=:time_sentinel_1_rtc, strata=strata)
+        fp_s2 = map(indices_s2) do band
+            arceme_bias_corrected_fp(band, ds, strata = strata)
+        end
+        fingerprint_sparse_s2 = setchunks(
+            YAXArrays.concatenatecubes(map(i -> i.fp, fp_s2), indexdim_s2),
+            (; band_s2 = length(indices_s2)),
+        )
+        fingerprint_uncor_sparse_s2 = setchunks(
+            YAXArrays.concatenatecubes(map(i -> i.fp_uncorrected, fp_s2), indexdim_s2),
+            (; band_s2 = length(indices_s2)),
+        )
+        #fingerprints_s2 = time_aggregate_fingerprint(fingerprint_sparse_s2, eventdate, indexdim_s2, :time_sentinel_2_l2a)
+        s2cubes = (
+            :s2_indices => fingerprint_sparse_s2,
+            :uncorrected_s2_indices => fingerprint_uncor_sparse_s2,
+            :class_fractions => fp_s2[1].class_fractions,
+        )
     end
-    fingerprint_sparse_s1 = setchunks(YAXArrays.concatenatecubes(map(i -> i.fp, fp_s1), indexdim_s1), (; band_s1=length(indices_s1)))
-
-    Dataset(; s2_indices=fingerprint_sparse_s2, s1_indices=fingerprint_sparse_s1, uncorrected_s2_indices=fingerprint_uncor_sparse_s2, class_fractions=fp_s1[1].class_fractions)
+    Dataset(; s1cubes..., s2cubes...)
 end
 
 function arceme_fractions(ds)
     c = counter(ds.ESA_LC[time=1].data[:, :])
     npix = length(ds.ESA_LC)
     #Find maximum landcover in both cubes
-    classkeys, classnames = collect(keys(arceme_legends["ESA_LC"])), collect(values(arceme_legends["ESA_LC"]))
-    classax = DD.Dim{:lc}(classnames)
+    classkeys, classnames =
+        collect(keys(arceme_legends["ESA_LC"])), collect(values(arceme_legends["ESA_LC"]))
+    classax = DD.Dim{:class}(classnames)
     fracs = zeros(Float64, length(classax))
     for (k, v) in c
         iclass = findfirst(==(k), classkeys)
         fracs[iclass] = v / npix
     end
-    ds.cubes[:lc_fraction] = YAXArray((classax,), fracs)
-    ds.axes[:lc] = classax
-    sumcl = xmap(ds.cloud_mask ⊘ (:x, :y), ds.SCL ⊘ (:x, :y), inplace=false, output=XOutput(outtype=Float64)) do cl, scl
+    ds.cubes[:class_fraction] = YAXArray((classax,), fracs)
+    ds.axes[:class] = classax
+    sumcl = xmap(
+        ds.cloud_mask ⊘ (:x, :y),
+        ds.SCL ⊘ (:x, :y),
+        inplace = false,
+        output = XOutput(outtype = Float64),
+    ) do cl, scl
         sum(i -> ARCEMEAnalysis._is_cloud(i...), zip(cl, scl))
     end
     cloudfrac = sumcl[1, 1, :].data ./ length(ds.x) ./ length(ds.y)
@@ -637,15 +886,74 @@ function arceme_fractions(ds)
     ds
 end
 
+"""
+Compute a cloud mask with a 5x5 window halo
+"""
+function arceme_cloudmask_halo(d)
+    combined_clouds = xmap(
+        d.cloud_mask,
+        d.SCL,
+        inplace = false,
+        output = XOutput(outtype = UInt8),
+    ) do cl, scl
+        ARCEMEAnalysis._is_clear_or_shadow(cl, scl) + 2 * (cl==3)
+    end
+    inwindows = (
+        DAE.InputArray(
+            combined_clouds.data,
+            windows = (
+                DAE.MovingWindow(-1, 1, 5, 1000, (1, 1000)),
+                DAE.MovingWindow(-1, 1, 5, 1000, (1, 1000)),
+                1:size(combined_clouds, 3),
+            ),
+        ),
+    )
+    outspec = (
+        DAE.create_outwindows(
+            (1000, 1000, size(combined_clouds, 3)),
+            chunks = (500, 500, 25),
+        ),
+    )
+    halofunc(x) = any(==(0), x) ? 0 : all(==(1), x) ? 1 : 2
+    op = DAE.GMDWop(
+        inwindows,
+        outspec,
+        DAE.create_userfunction(halofunc, UInt8, is_mutating = false),
+    )
+    d.cubes[:cloudmask_halo] = YAXArray(
+        d.cloud_mask.axes,
+        DAE.results_as_diskarrays(op)[1],
+        Dict("legend"=>"0 => Clouds or missing, 1=>clear sky, 2=> clear sky or shadow"),
+    )
+end
+
+function arceme_precompute_fingerprints(
+    tmpdir = "./fingerprints";
+    batch = "ARCEME-DC-DHP-GLOBAL",
+)
+    @showprogress pmap(Iterators.product(1:2, arceme_validpairs(; batch))) do (iev, ev_p)
+        ev = ev_p[iev]
+        ds = arceme_open(ev; batch)
+        fp = ARCEMEAnalysis.arceme_index_fingerprints(ds)
+        savedataset(fp, path = joinpath(tmpdir, arceme_cubename(ev)), overwrite = true)
+    end
+
+    outputzipfile = joinpath(ARCEMEAnalysis.local_cubepath, "$batch-fingerprints.zip")
+    isfile(outputzipfile) && rm(outputzipfile)
+    ARCEMEAnalysis.zip_dir(outputzipfile, tmpdir)
+end
+
 using Reexport: @reexport
 
 include("download.jl")
 include("spatialdebias.jl")
 include("spectral_helpers.jl")
 include("s1_helpers.jl")
-include("hrl_helpers.jl")
+include("hrl/hrl_helpers.jl")
 
 @reexport using .HRL
 
 include("plots/fingerprint_plots.jl")
+include("plots/cube_plots.jl")
+
 end #module
\ No newline at end of file
diff --git a/src/download.jl b/src/download.jl
index b364e1f..a32af0a 100644
--- a/src/download.jl
+++ b/src/download.jl
@@ -32,7 +32,7 @@ end
 
 const local_cubepath = @load_preference("arceme_localpath")
 
-function arceme_download_batch(batch="ARCEME-DC-6")
+function arceme_download_batch(batch="ARCEME-DC-DHP-GLOBAL"; subset=nothing)
     if local_cubepath === nothing
         error("You need to set the local cube path first. Please run `arceme_localpath(path)` first.")
     end
@@ -44,11 +44,16 @@ function arceme_download_batch(batch="ARCEME-DC-6")
     aresp = HTTP.get("$httpstore/$batch/",query=Dict("format"=>"json","delimiter"=>"/"))
     allarrays = map(i->strip(i["subdir"],'/'),JSON.parse(aresp.body))
 
+    if subset !== nothing
+        filter!(in(subset), allarrays)
+    end
+
     @showprogress for current_cube in allarrays
         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"))
+        isfile(outfilename) && rm(outfilename)
         ZipWriter(outfilename) do w
             for f in files_in_cube
                 resp = HTTP.get("$httpstore/$batch/$f");
diff --git a/src/hrl_helpers.jl b/src/hrl/hrl_helpers.jl
similarity index 64%
rename from src/hrl_helpers.jl
rename to src/hrl/hrl_helpers.jl
index aadd907..8cf6127 100644
--- a/src/hrl_helpers.jl
+++ b/src/hrl/hrl_helpers.jl
@@ -1,7 +1,7 @@
 module HRL
 using ARCEMEAnalysis
 import CondaPkg
-using  PythonCall: pyimport, pyconvert, pylist
+using PythonCall: pyimport, pyconvert, pylist
 import GeoJSON
 using Preferences: @set_preferences!, @load_preference
 import ArchGDAL as AG
@@ -26,7 +26,7 @@ end
 const hrl_localpath = @load_preference("hrl_localpath", pwd())
 
 if !isdir(hrl_localpath)
-  mkdir(hrl_localpath)
+    mkdir(hrl_localpath)
 end
 
 
@@ -37,18 +37,20 @@ Sets up client to Harmonised Data Access API on WEkEO
 Returns HDA client as Python Object
 """
 function hrl_hda()
-# Harmonised Data Access API on WEkEO
-  hda = try pyimport("hda") #("Client", "Configuration")
-  catch
-      CondaPkg.add("hda")
-      pyimport("hda")
-  end
-  # Configure user's credentials with a .hdarc file
-  hda_client = try hda.Client()
-  catch
-      @error "Please configure your WEkEO user's credentials with a .hdarc file  (see https://hda.readthedocs.io/en/latest/usage.html#client-configuration for details)"
-  end
-  return hda_client
+    # Harmonised Data Access API on WEkEO
+    hda = try
+        pyimport("hda") #("Client", "Configuration")
+    catch
+        CondaPkg.add("hda")
+        pyimport("hda")
+    end
+    # Configure user's credentials with a .hdarc file
+    hda_client = try
+        hda.Client()
+    catch
+        @error "Please configure your WEkEO user's credentials with a .hdarc file  (see https://hda.readthedocs.io/en/latest/usage.html#client-configuration for details)"
+    end
+    return hda_client
 end
 
 """
@@ -58,12 +60,12 @@ Search for all datasets available on WEkEO in whose name `filteron` occurs. Very
 
 Returns a `Vector` of `String`s
 """
-function hda_datasets(;filteron="HRL")
-  hda_client = hrl_hda()
-  datasets = hda_client.datasets() 
-  dataset_ids = [item["dataset_id"] for item in datasets] 
-  # convert to julia
-  pyconvert(Vector{String}, dataset_ids) |> filter(x -> occursin(filteron, x))
+function hda_datasets(; filteron = "HRL")
+    hda_client = hrl_hda()
+    datasets = hda_client.datasets()
+    dataset_ids = [item["dataset_id"] for item in datasets]
+    # convert to julia
+    pyconvert(Vector{String}, dataset_ids) |> filter(x -> occursin(filteron, x))
 end
 
 """
@@ -73,9 +75,12 @@ Lists all valid productType for `dataset_id` on WEkEO.
 Returns a vector of strings
 """
 function hda_products(dataset_id)
-  hda_client = hrl_hda()
-  metadata = hda_client.metadata(dataset_id)
-  map(x -> x["const"], pyconvert(Vector{Dict},metadata["properties"]["productType"]["oneOf"]))
+    hda_client = hrl_hda()
+    metadata = hda_client.metadata(dataset_id)
+    map(
+        x -> x["const"],
+        pyconvert(Vector{Dict}, metadata["properties"]["productType"]["oneOf"]),
+    )
 end
 
 """
@@ -84,10 +89,10 @@ end
 Lists all valid years or epochs for `dataset_id` on WEkEO.
 Returns a vector of strings
 """
-function  hda_years(dataset_id)
-  hda_client = hrl_hda()
-  metadata = hda_client.metadata(dataset_id)
-  map(x -> x["const"], pyconvert(Vector{Dict},metadata["properties"]["year"]["oneOf"]))
+function hda_years(dataset_id)
+    hda_client = hrl_hda()
+    metadata = hda_client.metadata(dataset_id)
+    map(x -> x["const"], pyconvert(Vector{Dict}, metadata["properties"]["year"]["oneOf"]))
 end
 
 
@@ -106,20 +111,26 @@ Returns a Vector of filenames.
 
 See also `hda_datasets`, `hda_products`, `hda_years`
 """
-function hrl_download(dataset_id::String, coord::Tuple, yr::Union{String,Int}; extent=0.05, product::Union{Nothing,String}=nothing)
+function hrl_download(
+    dataset_id::String,
+    coord::Tuple,
+    yr::Union{String,Int};
+    extent = 0.05,
+    product::Union{Nothing,String} = nothing,
+)
     hda_client = hrl_hda()
     json = pyimport("json")
 
-    (lon,lat) = coord
+    (lon, lat) = coord
     fnames = []
 
     # Build query as a JSON object
     if !isnothing(product)
-      pt = """,
-      "productType": "$(product)"
-      """
+        pt = """,
+        "productType": "$(product)"
+        """
     else
-      pt = ""
+        pt = ""
     end
 
     query = """
@@ -137,35 +148,43 @@ function hrl_download(dataset_id::String, coord::Tuple, yr::Union{String,Int}; e
       "startIndex": 0$(pt)
     }
     """
-    
+
     matches = hda_client.search(json.loads(query))
     try
-    for item in matches.results
-      # println(item["id"])
-      filename = pyconvert(String, item["id"])*".zip"
-      # store info:
-      push!(fnames, filename)
-    end
+        for item in matches.results
+            # println(item["id"])
+            filename = pyconvert(String, item["id"]) * ".zip"
+            # store info:
+            push!(fnames, filename)
+        end
     catch
-      @warn "No matches found for $(dataset_id) at $(coordinates) in year $(yr)"
-      return fnames
+        @warn "No matches found for $(dataset_id) at $(coordinates) in year $(yr)"
+        return fnames
     end
 
     # Download
-    todownload = [!isfile(joinpath(hrl_localpath, pyconvert(String, item["id"])*".zip")) for  item in matches.results]
-    @showprogress desc="Downloading..." for i in pylist(findall(todownload).-1)
-    #   println("Downloading $(matches.results[i]["id"])")
-      matches[i].download(download_dir=hrl_localpath) 
+    todownload = [
+        !isfile(joinpath(hrl_localpath, pyconvert(String, item["id"]) * ".zip")) for
+        item in matches.results
+    ]
+    @showprogress desc = "Downloading..." for i in pylist(findall(todownload) .- 1)
+        #   println("Downloading $(matches.results[i]["id"])")
+        matches[i].download(download_dir = hrl_localpath)
     end
     return fnames
 end
 
-function hrl_download(datasets::Vector{String}, coord::Tuple, yr::Union{String,Int}; kwargs...)
-  fnames=[];
-  for dataset_id in datasets
-      append!(fnames, HRL.hrl_download(dataset_id, coord, yr; kwargs...))
-  end
-  return fnames
+function hrl_download(
+    datasets::Vector{String},
+    coord::Tuple,
+    yr::Union{String,Int};
+    kwargs...,
+)
+    fnames = []
+    for dataset_id in datasets
+        append!(fnames, HRL.hrl_download(dataset_id, coord, yr; kwargs...))
+    end
+    return fnames
 end
 
 const hrl_dict = Dict(
@@ -200,7 +219,7 @@ const hrl_dict = Dict(
     "GRACL" => "Grassland Confidence Layer",
     "HER" => "Herbaceous Cover",
     "PLOUGH" => "Ploughing Indicator",
-    "GRAME"  => "Grassland Mowing Events",
+    "GRAME" => "Grassland Mowing Events",
     "GRAMECL" => "Grassland Mowing Events Confidence Layer",
     "GRAMD1" => "Grassland Mowing Dates of first mowing event",
     "GRAMD2" => "Grassland Mowing Dates of second mowing event",
@@ -209,13 +228,9 @@ const hrl_dict = Dict(
 )
 
 const hrl_legends = Dict(
-    "MCTY" => SortedDict(
-        1 => "No cropland",
-        2 => "Annual crop",
-        3 => "Permanent crop"
-    ),
+    "MCTY" => SortedDict(1 => "No cropland", 2 => "Annual crop", 3 => "Permanent crop"),
     "CTY" => SortedDict(
-        0 =>  "No cropland",
+        0 => "No cropland",
         1110 => "Wheat",
         1120 => "Barley",
         1130 => "Maize",
@@ -235,7 +250,7 @@ const hrl_legends = Dict(
         2320 => "Nuts",
         3100 => "Unclassified annual crop",
         3200 => "Unclassified permanent crop",
-        65535 => "outside area"
+        65535 => "outside area",
     ),
     "CTYCL" => "0-100 Probability expressed as a percentage. 253 No Cropland. 255 Outside area",
     "CPMCH" => "Harvest date as YYDOY where YY = last 2 digits of the year (e.g., 19 for 2019) and DOY is the day of the year (1-366). 
@@ -254,7 +269,7 @@ const hrl_legends = Dict(
         65532 No cropping season detected. 
         65533 Growing season extends beyond timeframe. 
         64534 No confidence could be calculated.
-        65535 Outside area.",   
+        65535 Outside area.",
     "CPMCE" => "Emergence date as YYDOY where YY = last 2 digits of the year (e.g., 19 for 2019) and DOY is the day of the year (1-366). 
         0 No annual cropland. 
         65526 Fallow land. 
@@ -328,18 +343,18 @@ const hrl_legends = Dict(
         64534 No confidence could be calculated.
         65535 Outside area.",
     "CPSCT" => SortedDict(
-        0 => "No annual cropland", 
+        0 => "No annual cropland",
         1 => "Short Summer",
         2 => "Long Summer",
         3 => "Short Winter",
-        4 => "Long Winter", 
-        65526 => "Fallow land", 
-        65527 => "No cropping pattern detected",  
+        4 => "Long Winter",
+        65526 => "Fallow land",
+        65527 => "No cropping pattern detected",
         65530 => "No secondary crop growing season delineated",
-        65531 => "Not enough data", 
-        65532 => "No cropping season detected",  
+        65531 => "Not enough data",
+        65532 => "No cropping season detected",
         65533 => "Growing season extends beyond timeframe",
-        65535 => "Outside area"
+        65535 => "Outside area",
     ),
     "CPSCE" => "Emergence date as YYDOY where YY = last 2 digits of the year (e.g., 19 for 2019) and DOY is the day of the year (1-366). 
         0 No annual cropland. 
@@ -369,22 +384,19 @@ const hrl_legends = Dict(
         65533 Growing season extends beyond timeframe. 
         64534 No confidence could be calculated.
         65535 Outside area.",
-    "CPFLP" => SortedDict(
-        0 => "No fallow land",
-        1 => "Fallow land",
-        65535 => "Outside area"
-    ),
+    "CPFLP" =>
+        SortedDict(0 => "No fallow land", 1 => "Fallow land", 65535 => "Outside area"),
     "CPFLPCL" => "0-100 Probability expressed as a percentage. 253 No Cropland. 65535 Outside area",
     "CPCSY" => SortedDict(
-        0 => "No annual cropland", 
+        0 => "No annual cropland",
         1 => "One growing season",
         2 => "Two growing seasons",
-        65526 => "Fallow land", 
-        65527 => "No cropping pattern detected",  
-        65531 => "Not enough data", 
-        65532 => "No cropping season detected",  
+        65526 => "Fallow land",
+        65527 => "No cropping pattern detected",
+        65531 => "Not enough data",
+        65532 => "No cropping season detected",
         65533 => "Growing season extends beyond timeframe",
-        65535 => "Outside area"
+        65535 => "Outside area",
     ),
     # Forest
     "DLT" => SortedDict(
@@ -405,7 +417,7 @@ const hrl_legends = Dict(
     # Grassland
     "GRA" => SortedDict(
         0 => "all non-grassland areas",
-        1 => "grassland", 
+        1 => "grassland",
         255 => "outside area",
     ),
     "GRACL" => "0-100: Classification confidence. 253: All non-grassland areas. 255: outside area",
@@ -419,14 +431,14 @@ const hrl_legends = Dict(
     100: Change in herbaceous cover. 
     253: no ploughing information. 
     255: outside area",
-    "GRAME"  => SortedDict(
-        0 => "no mowing detected", 
+    "GRAME" => SortedDict(
+        0 => "no mowing detected",
         1 => "1 mowing event detected",
         2 => "2 mowing events",
         3 => "3 mowing events",
         4 => "4 mowing events",
         253 => "all non-herbaceous areas",
-        255 => "outside area"
+        255 => "outside area",
     ),
     "GRAMECL" => "0-100: Mowing detection confidence. 253: All non-mowing areas. 255: outside area",
     "GRAMD1" => "0: no mowing detected. 1-366: Start (DOY) of each mowing event. 65533: all non-herbaceous areas. 65535: outside area",
@@ -443,111 +455,145 @@ const hrl_legends = Dict(
 reprojects and clips files in `fnames` to match `event`'s coordinate reference system and extent,
  and saves them as a Zarr group in `batch-HRL` under `ARCEMEAnalysis.local_cubepath`.
 """
-function hrl_warp(cubename, fnames::Vector; batch="ARCEME-DC-6")
-  ds = arceme_open(cubename; batch)
-  event_date = event_date = arceme_eventdate(ds)
-  tmp = tempname(;suffix=".tif")
-  res = ds.properties["resolution"]
-  tmp = tempname(;suffix=".tif")
-
-  for hrl in keys(hrl_dict)
-    hrl_fnames = fnames[map(contains("$(hrl)_"), fnames)]
-    if isempty(hrl_fnames)
-        continue
+function hrl_warp(cubename, fnames::Vector; batch = "ARCEME-DC-6", suffix = "", zip = true)
+    ds = arceme_open(cubename; batch)
+    event_date = arceme_eventdate(ds)
+    res = ds.properties["resolution"]
+    tmp = tempname(; suffix = ".tif")
+
+    for hrl in keys(hrl_dict)
+        hrl_fnames = fnames[map(contains("$(hrl)_"), fnames)]
+        if isempty(hrl_fnames)
+            continue
+        end
+        run(
+            Cmd(
+                `gdalwarp 
+                    -t_srs $(AG.toPROJ4(AG.importEPSG(ds.properties["epsg"]))) 
+                    -tr $res $res
+                    -te $(minimum(ds.x)) $(minimum(ds.y)) $(maximum(ds.x)+res) $(maximum(ds.y)+res) 
+                    -r near
+                    $(["/vsizip/$(joinpath(hrl_localpath, hrl_fnames[i], "$(split(hrl_fnames[i],".")[1]).tif") )" for i in eachindex(hrl_fnames)]) 
+                    $tmp 
+                    -overwrite 
+                    `,
+            ),
+        )
+        c = Cube(tmp)
+        # write to zarr
+        props = Dict(
+            "long_name" => hrl_dict["$(split(hrl_fnames[1],"_")[3])"],
+            "year_coverage" => split(hrl_fnames[1], "_")[4][2:5], #year(event_date),
+            "legend" => hrl_legends["$(split(hrl_fnames[1],"_")[3])"],
+            "epsg" => ds.properties["epsg"],
+            "resolution" => 10,
+            "unit" => "meter",
+            "source" => "$(join(hrl_fnames, ", ")) © European Union, Copernicus Land Monitoring Service 2025, European Environment Agency (EEA)",
+            "processing" => "Downloaded from wEkEO. Reprojected to EPSG $(ds.properties["epsg"]) using nearest neighour, mosaicked and clipped to AOI with gdalwarp",
+            "processed_by" => "M. Weynants, MPI-BGC",
+            "processed_date" => string(Date(now())),
+        )
+        c2 = YAXArray((ds.x, ds.y), c.data, props)
+        nt = (Symbol("$(split(hrl_fnames[1],"_")[3])$suffix") => c2,)
+        tds = Dataset(; nt...)
+        tds = setchunks(tds, (500, 500))
+        savedataset(
+            tds,
+            path = joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", cubename),
+            append = true,
+        )
     end
-    run(Cmd(
-        `gdalwarp 
-            -t_srs $(AG.toPROJ4(AG.importEPSG(ds.properties["epsg"]))) 
-            -tr $res $res
-            -te $(minimum(ds.x)) $(minimum(ds.y)) $(maximum(ds.x)+res) $(maximum(ds.y)+res) 
-            -r near
-            $(["/vsizip/$(joinpath(hrl_localpath, hrl_fnames[i], "$(split(hrl_fnames[i],".")[1]).tif") )" for i in eachindex(hrl_fnames)]) 
-            $tmp 
-            -overwrite 
-            `))
-    c = Cube(tmp)
-    # write to zarr
-    props = Dict(
-        "long_name" => hrl_dict["$(split(hrl_fnames[1],"_")[3])"],
-        "year_coverage" => year(event_date),
-        "legend" => hrl_legends["$(split(hrl_fnames[1],"_")[3])"],
-        "epsg" => ds.properties["epsg"], 
-        "resolution" => 10,
-        "unit" => "meter",
-        "source" => "$(join(hrl_fnames, ", ")) © European Union, Copernicus Land Monitoring Service 2025, European Environment Agency (EEA)",
-        "processing" => "Downloaded from wEkEO. Reprojected to EPSG $(ds.properties["epsg"]) using nearest neighour, mosaicked and clipped to AOI with gdalwarp",
-        "processed_by" => "M. Weynants, MPI-BGC",
-        "processed_date" => string(Date(now())),
+    # zip ?
+    !zip && return joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", cubename)
+
+    isfile(
+        joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", string(cubename, ".zip")),
+    ) && rm(
+        joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", string(cubename, ".zip")),
     )
-    c2 = YAXArray((ds.x,ds.y), c.data, props, )
-    nt = (Symbol("$(split(hrl_fnames[1],"_")[3])")=>c2,)
-    tds = Dataset(;nt...)
-    tds = setchunks(tds, (500,500))
-    savedataset(tds, 
-        path=joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", "$(split(cubename,".")[1]).zarr"), 
-        append=true)
-    end
-    # zip
-    run(Cmd(`zip -0 -r ../$(split(cubename,".")[1]).zarr.zip .`, dir=joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", "$(split(cubename,".")[1]).zarr")))
-    rm(joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", "$(split(cubename,".")[1]).zarr"), recursive=true)
-    return joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", "$(split(cubename,".")[1]).zarr.zip")
+    ARCEMEAnalysis.zip_dir(
+        joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", "$cubename.zip"),
+        joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", cubename),
+    )
+    rm(joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", cubename), recursive = true)
+    return joinpath(ARCEMEAnalysis.local_cubepath, "$(batch)-HRL", "$cubename.zip")
 end
 
-hrl_warp(event::ARCEMEAnalysis.Event; kwargs...) = hrl_warp(arceme_cubename(event); kwargs...)
-
-function hrl_warp(cubename::String; batch="ARCEME-DC-6", dataset_id::Union{String, Vector{String}}=["EO:EEA:DAT:HRL:CRL", "EO:EEA:DAT:HRL:TCF", "EO:EEA:DAT:HRL:GRA"])
-  ds = arceme_open(cubename; batch=batch)
-  coord = arceme_coordinates(ds)
-  yr = year(arceme_eventdate(ds))
-  fnames = hrl_download(dataset_id, coord, yr)
-  hrl_warp(cubename, fnames; batch)
+hrl_warp(event::ARCEMEAnalysis.Event; kwargs...) =
+    hrl_warp(arceme_cubename(event); kwargs...)
+
+function hrl_warp(
+    cubename::String;
+    batch = "ARCEME-DC-6",
+    dataset_id::Union{String,Vector{String}} = [
+        "EO:EEA:DAT:HRL:CRL",
+        "EO:EEA:DAT:HRL:TCF",
+        "EO:EEA:DAT:HRL:GRA",
+    ],
+    yr = nothing,
+    kwargs...,
+)
+    ds = arceme_open(cubename; batch)
+    coord = arceme_coordinates(ds)
+    if isnothing(yr)
+        yr = year(arceme_eventdate(ds))
+    end
+    fnames = hrl_download(dataset_id, coord, yr)
+    hrl_warp(cubename, fnames; batch, kwargs...)
 end
 
 """
      arceme_stats(ds::YAXArrays.Dataset, name; classes=nothing)
-     arceme_stats(ev::Event, name; batch="ARCEME-DC-6", classes=nothing)
+     arceme_stats(ev::Event, name; batch="ARCEME-DC-8", classes=nothing)
 
 Get the  statistics of layer `name` for the ARCEME data cube `ds`.
 Default classes are extracted from `hrl_legends[name]`. Alternative classes can be provided as a dictionary.
 """
-function arceme_stats(ds, name; classes=nothing)
+function arceme_stats(ds, name; classes = nothing)
     if isnothing(classes)
         classes = hrl_legends[name]
     end
-    cdr = counter(ds[name]);
-    map(collect(classes)) do (k,v)
-        count=get(cdr, k, 0)
-        (key=k, class=v, count=count, fraction=count/1000000)
+    cdr = counter(ds[name])
+    map(collect(classes)) do (k, v)
+        count = get(cdr, k, 0)
+        (
+            key = k,
+            class = v,
+            count = count,
+            fraction = count / (length(ds.x) * length(ds.y)),
+        )
     end
 end
-arceme_stats(ev::ARCEMEAnalysis.Event, name; batch="ARCEME-DC-6", classes=nothing) = arceme_stats(arceme_open(ev, batch=batch), name; classes = classes)
+arceme_stats(ev::ARCEMEAnalysis.Event, name; batch = "ARCEME-DC-8", classes = nothing) =
+    arceme_stats(arceme_open(ev, batch = batch), name; classes = classes)
 
 function ctykeymap(k)
     k == 0 && return 1
     k > 3000 && return 1
     ctL1 = k ÷ 1000
-    ctL2 = (k - 1000*ctL1) ÷ 100
-    ctL3 = (k - 1000*ctL1 - 100*ctL2) ÷ 10
-    ctL1==1 && ctL2==1 && return ctL3+1 # cereals
-    ctL1==1 && ctL2==2 && return ctL3+6 # vegetables and pulses
-    ctL1==1 && ctL2==3 && return ctL3+8 # roots (potato, sugar beet)
-    ctL1==1 && ctL2==4 && return ctL3+10 # oil (sunflower, soy, rapeseed)
-    ctL1==2 && ctL2==1 && return ctL3+15 # fiber (flax, cotton, hemp)
-    ctL1==2 && ctL2==2 && return ctL2+14 # grape, olive
-    ctL1==2 && ctL2==3 && return ctL3+16 # fruit, nuts
+    ctL2 = (k - 1000 * ctL1) ÷ 100
+    ctL3 = (k - 1000 * ctL1 - 100 * ctL2) ÷ 10
+    ctL1 == 1 && ctL2 == 1 && return ctL3 + 1 # cereals
+    ctL1 == 1 && ctL2 == 2 && return ctL3 + 6 # vegetables and pulses
+    ctL1 == 1 && ctL2 == 3 && return ctL3 + 8 # roots (potato, sugar beet)
+    ctL1 == 1 && ctL2 == 4 && return ctL3 + 10 # oil (sunflower, soy, rapeseed)
+    ctL1 == 2 && ctL2 == 1 && return ctL3 + 15 # fiber (flax, cotton, hemp)
+    ctL1 == 2 && ctL2 == 2 && return ctL2 + 14 # grape, olive
+    ctL1 == 2 && ctL2 == 3 && return ctL3 + 16 # fruit, nuts
 end
 
 function mctykeymap(k)
     k == 0 && return 1 # non crop
     k == 65535 && return 1
     ctL1 = k ÷ 1000
-    ctL2 = (k - 1000*ctL1) ÷ 100
-    ctL3 = (k - 1000*ctL1 - 100*ctL2) ÷ 10
+    ctL2 = (k - 1000 * ctL1) ÷ 100
+    ctL3 = (k - 1000 * ctL1 - 100 * ctL2) ÷ 10
     ctL1 == 1 && return 2 # annual crops
     # ctL1 == 2 && ctL2 == 1 && return 2 # annual fiber
     ctL1 == 2 && return 3 # permanent crop
     k == 3100 && return 2 # Unclassified annual crop
     k == 3200 && return 3 # Unclassified permanent crop
 end
-export hrl_download, hrl_set_localpath, hda_datasets, hda_products, hda_years, hrl_warp, arceme_stats
+export hrl_download,
+    hrl_set_localpath, hda_datasets, hda_products, hda_years, hrl_warp, arceme_stats
 end #module
\ No newline at end of file
diff --git a/src/plots/cube_plots.jl b/src/plots/cube_plots.jl
new file mode 100644
index 0000000..a87b8ea
--- /dev/null
+++ b/src/plots/cube_plots.jl
@@ -0,0 +1,882 @@
+@reexport module CubePlots
+using Colors: Colorant
+using Makie
+using YAXArrays
+using Dates, DateFormats
+using ..ARCEMEAnalysis
+using ..ARCEMEAnalysis: local_cubepath, arceme_legends, Event, arceme_open
+using DataStructures: counter, OrderedDict
+using Statistics: mean, median
+using GeoMakie
+import DimensionalData as DD
+using Random: seed!
+using DataFrames
+using Interpolations
+using HypothesisTests: OneWayANOVATest
+using RollingFunctions
+
+
+
+ESA_LC_colormap = [
+    (color = "#006400", alpha = 255, value = 10, label = "Tree cover"),
+    (color = "#ffbb22", alpha = 255, value = 20, label = "Shrubland"),
+    (color = "#ffff4c", alpha = 255, value = 30, label = "Grassland"),
+    (color = "#f096ff", alpha = 255, value = 40, label = "Cropland"),
+    (color = "#fa0000", alpha = 255, value = 50, label = "Built-up"),
+    (color = "#b4b4b4", alpha = 255, value = 60, label = "Bare / sparse vegetation"),
+    (color = "#f0f0f0", alpha = 255, value = 70, label = "Snow and ice"),
+    (color = "#0064c8", alpha = 255, value = 80, label = "Permanent water bodies"),
+    (color = "#0096a0", alpha = 255, value = 90, label = "Herbaceous wetland"),
+    (color = "#00cf75", alpha = 255, value = 95, label = "Mangroves"),
+    (color = "#fae6a0", alpha = 255, value = 100, label = "Moss and lichen"),
+]
+clc_clrs = [
+    colorant"#006400",
+    colorant"#ffbb22",
+    colorant"#ffff4c",
+    colorant"#f096ff",
+    colorant"#fa0000",
+    colorant"#b4b4b4",
+    colorant"#f0f0f0",
+    colorant"#0064c8",
+    colorant"#0096a0",
+    colorant"#00cf75",
+    colorant"#fae6a0",
+]
+
+cty_clrs = cgrad(
+    vcat(
+        cgrad(:lipariS, 14, categorical = true).colors.colors,
+        cgrad(:bamakoS, 4, categorical = true).colors.colors,
+    ),
+    18,
+    categorical = true,
+)
+
+ctykeys = hcat(
+    values(arceme_legends["CTY"]) |> collect,
+    keys(arceme_legends["CTY"]) |> collect,
+    ARCEMEAnalysis.HRL.mctykeymap.(keys(arceme_legends["CTY"]) |> collect),
+    ARCEMEAnalysis.HRL.ctykeymap.(keys(arceme_legends["CTY"]) |> collect),
+)
+
+function getticks(mn, mx)
+    n = length(range(mn, mx))
+    collect((mn+(mx-mn)/(2*n)):((mx-mn)/n):mx)
+end
+
+function plotlc(
+    current_event::Event;
+    batch = "ARCEME-DC-8",
+    savef = true,
+    fname = "fig_ESA_LC_$(current_event.uid).png",
+)
+    ds = arceme_open(current_event, batch = batch)
+    f, ax, p = image(
+        ds.ESA_LC[time=1],
+        colormap = cgrad(clc_clrs; categorical = true),
+        colorrange = (10, 110),
+        # label = "Land cover",
+        axis = (title = "$(current_event.uid) - Land Cover", aspect = DataAspect()),
+    )
+    hidedecorations!(ax)
+    Colorbar(
+        f[1, 2],
+        colormap = cgrad(clc_clrs; categorical = true),
+        colorrange = (1, 11),
+        ticks = (getticks(1, 11), collect(values(arceme_legends["ESA_LC"]))[2:end]),
+    )
+    savef && save(fname, f)
+    f
+end
+
+function plotcty(
+    current_event::Event;
+    batch = "ARCEME-DC-8",
+    savef = true,
+    fname = "fig_CTY_$(current_event.uid).png",
+)
+    ds = arceme_open(current_event, batch = batch)
+    cty = map(ARCEMEAnalysis.HRL.ctykeymap, ds.CTY)
+    f, ax, p = image(
+        cty,
+        colormap = cgrad(cty_clrs; categorical = true),
+        colorrange = (1, 18),
+        # lowclip = :white,
+        highclip = :transparent,
+        axis = (title = "$(current_event.uid) - Main Crop Type", aspect = DataAspect()),
+    )
+    hidedecorations!(ax)
+    Colorbar(
+        f[1, 2],
+        colormap = cgrad(cty_clrs; categorical = true),
+        colorrange = (1, 18),
+        # lowclip = :white,
+        # highclip = :white,
+        ticks = (getticks(1, 18), [k for k in ctykeys[1:18, 1]]),
+    )
+    savef && save(fname, f)
+    f
+end
+
+yydoy2date(x::Integer) = DateTime("$(Int(2000 + x÷1e3))", "YY") + Day(x % 1e3)
+yydoy2date(x::Missing) = missing
+yydoy2date(x::Vector) = broadcast(yydoy2date, x)
+
+yydoy2ydec(x::Integer) = yeardecimal(yydoy2date(x))
+yydoy2ydec(x::Missing) = missing
+yydoy2ydec(x::Vector) = broadcast(yydoy2ydec, x)
+
+rep(A) = replace(
+    convert(Vector{Union{Missing,Int32}}, A),
+    0 => missing,
+    65526 => missing,
+    65527 => missing,
+    65530 => missing,
+    65531 => missing,
+    65532 => missing,
+    65533 => missing,
+    65534 => missing,
+    65535 => missing,
+)
+
+
+function plotmceh(
+    current_event::Event;
+    batch = "ARCEME-DC-8",
+    savef = true,
+    fname = "fig_CPMCEH_$(current_event.uid).png",
+)
+    ds = arceme_open(current_event, batch = batch)
+    event_date = arceme_eventdate(current_event)
+    mce = ds.CPMCE.data[:] |> rep |> yydoy2ydec
+    mch = ds.CPMCH.data[:] |> rep |> yydoy2ydec
+
+    fig = Figure(size = (1000, 500))
+    ax = Axis(
+        fig[1, 1],
+        title = current_event.uid,
+        xlabel = "left: CPMCE / right: CPMCH",
+        ylabel = "Year decimal",
+    )
+    cty = ARCEMEAnalysis.HRL.ctykeymap.(ds.CTY.data[:])[.!(ismissing.(mce))]
+    for i = 1:14
+        hist!(
+            ax,
+            mce[.!(ismissing.(mce))][cty .== i],
+            scale_to = -0.5,
+            offset = i,
+            direction = :x,
+            color = cty_clrs.colors.colors[i],
+            label = "$i: $(ctykeys[i,1])",
+        )
+        hist!(
+            ax,
+            mch[.!(ismissing.(mch))][cty .== i],
+            scale_to = 0.5,
+            offset = i,
+            direction = :x,
+            color = cty_clrs.colors.colors[i],
+        )
+    end
+    ax.xticks = 1:14
+    # Colorbar(fig[1,2], label = "Main crop type",
+    #     colormap = cgrad(:lipariS, 14, categorical = true),
+    #     colorrange = (1,14),
+    #     ticks=(getticks(1,14),ctykeys[1:14,1]),
+    # )
+    hlines!(
+        ax,
+        yeardecimal(event_date),
+        color = "grey",
+        label = "Event date: $(Date(event_date))",
+    )
+    fig[1, 2] = Legend(fig, ax, "Main Crop Types", framevisible = false)
+
+    savef && save(fname, fig)
+    fig
+end
+
+function piecty(
+    current_event::Event;
+    batch = "ARCEME-DC-8",
+    savef = true,
+    fname = "fig_CTY_pie_$(current_event.uid).png",
+)
+    cubename = arceme_cubename(current_event)
+    fp_cty = open_dataset(
+        joinpath(
+            ARCEMEAnalysis.local_cubepath,
+            "$(batch)-fingerprints",
+            "$(split(cubename,".")[1])__fp_CTY.zarr",
+        ),
+    )
+    f, ax, plt = pie(
+        fp_cty.class_fractions.data[2:18],
+        color = cty_clrs[2:end],
+        radius = 4,
+        inner_radius = 2,
+        strokecolor = :white,
+        strokewidth = 2,
+        axis = (autolimitaspect = 1,),
+        label = [k => (; color = c) for (k, c) in zip(ctykeys[2:18, 1], cty_clrs[2:end])],
+    )
+    leg = Legend(f[1, 2], ax)
+    hidedecorations!(ax)
+
+    savef && save(fname, f)
+    f
+end
+
+function fpfig(
+    current_event::Event,
+    class;
+    idx = "kNDVI",
+    mclass = "Annual crop",
+    batch = "ARCEME-DC-8",
+    savef = true,
+    fname = "fig_fp_$(class)_$(current_event.uid).png",
+)
+    cubename = arceme_cubename(current_event)
+    event_date = arceme_eventdate(current_event)
+    fp_cty = open_dataset(
+        joinpath(
+            ARCEMEAnalysis.local_cubepath,
+            "$(batch)-fingerprints",
+            "$(split(cubename,".")[1])__fp_CTY.zarr",
+        ),
+    )
+
+    fig = Figure()
+    tempo_s2 = lookup(fp_cty.time_sentinel_2_l2a)
+    tick_positions = datetime2unix.((event_date-Year(1)):Month(4):(event_date+Year(1)))  # Convert dates to Unix time (positions)
+    tick_labels = [
+        "$(Dates.format(t, "yyyy-mm-dd") )" for
+        t = (event_date-Year(1)):Month(4):(event_date+Year(1))
+    ]
+    n = length(tempo_s2)
+    ax = Axis(
+        fig[1, 1],
+        title = "$(current_event.uid)",
+        xlabel = "Time",
+        ylabel = idx,
+        xticks = (tick_positions, tick_labels),
+    )
+    if !isnothing(mclass)
+        mfp_cty = open_dataset(
+            joinpath(
+                ARCEMEAnalysis.local_cubepath,
+                "$(batch)-fingerprints",
+                "$(split(cubename,".")[1])__fp_MCTY.zarr",
+            ),
+        )
+        scatterlines!(
+            ax,
+            datetime2unix.(tempo_s2),
+            mfp_cty.uncorrected_s2_indices[band_s2=At(idx), class=At(mclass)].data[:],
+            label = "$(mclass) - $idx",
+            linewidth = 2,
+        )
+    end
+    scatterlines!(
+        ax,
+        datetime2unix.(tempo_s2),
+        fp_cty.uncorrected_s2_indices[band_s2=At(idx), class=At(class)].data[:],
+        label = "$(class) - $idx",
+        linewidth = 1,
+    )
+    # ax.xticks = ([1,n÷2+1,n], string.(Date.(tempo_s2[[1,n÷2+1,n]])))
+    fig[2, 1] = Legend(fig, ax, framevisible = false, orientation = :horizontal, nbanks = 1)
+
+    savef && save(fname, fig)
+    fig
+end
+
+function getsplines(current_event; bands = ["kNDVI", "NDVI"], batch = "ARCEME-DC-8")
+    ds = arceme_open(current_event, batch = batch)
+    firstday = arceme_eventdate(current_event) - Year(1)
+    lastday = arceme_eventdate(current_event) + Year(1)
+    nms = Millisecond(lastday - firstday).value
+    linscal = d -> Millisecond(d - firstday).value * 2 / nms
+    x_s2 = map(linscal, ds.time_sentinel_2_l2a.val)
+    outax = YAXArrays.time(range(firstday, lastday, step = Day(1)))
+    s2 = YAXArray(
+        (ds.x, ds.y, ds.time_sentinel_2_l2a, Dim{:band_s2}(bands)),
+        cat(tuple([ds[bands[i]].data for i in eachindex(bands)]...)..., dims = 4),
+    )
+    r_s2 = xmap(
+        ARCEMEAnalysis.FingerprintPlots.filter_spline,
+        s2 ⊘ :time_sentinel_2_l2a,
+        output = XOutput(outax),
+        function_args = (x_s2, length(outax)),
+    )
+    return r_s2
+end
+
+function plot_cpbsb_splines(
+    current_event::Event,
+    crop;
+    batch = "ARCEME-DC-8",
+    idx = "kNDVI",
+    nsamples = 100,
+    nmonths = 6,
+    random = true,
+    ind = [1, 10, 100, 1000],
+    savef = true,
+    fname = "fig_$(idx)_splines_$(crop)_cpbsb_$(current_event.uid).png",
+)
+    ds = arceme_open(current_event, batch = batch)
+    r_s2 = getsplines(current_event)
+    event_date = arceme_eventdate(current_event)
+
+    cls = ctykeys[findfirst(ctykeys[:, 1] .== crop), 2]
+    indcrop = findall(x -> x == cls, ds.CTY.data[:, :])
+    f = Figure()
+    tick_positions =
+        datetime2unix.((event_date-Year(1)):Month(nmonths):(event_date+Year(1)))  # Convert dates to Unix time (positions)
+    tick_labels = [
+        "$(Dates.format(t, "yyyy-mm-dd") )" for
+        t = (event_date-Year(1)):Month(nmonths):(event_date+Year(1))
+    ]
+    a1 = Axis(
+        f[1, 1],
+        title = "$(current_event.uid) - $crop",
+        xlabel = "Date",
+        ylabel = idx,
+        xticks = (tick_positions, tick_labels),
+    )
+    if random
+        seed!(42)
+        cind = rand(indcrop, nsamples)
+    else
+        cind = indcrop[ind]
+    end
+    for i in cind
+        cpbsb = ds.CPBSB.data[i[1], i[2]]
+        cpbsb > 65000 && continue
+        lines!(
+            a1,
+            datetime2unix.(lookup(r_s2.time)),
+            r_s2[band_s2=At(idx)].data[:, i[1], i[2], 1],
+            color = cpbsb,
+            colorrange = (0, 295),
+            colormap = (:viridis, 0.5),
+        )
+    end
+    Colorbar(
+        f[1, 2],
+        colorrange = (0, 295),
+        colormap = (:viridis, 0.5),
+        label = "Bare soil before crop emergence (days)",
+    )
+    savef && save(fname, f)
+    f
+end
+
+function getindcrop(current_event, crop; batch = "ARCEME-DC-8")
+    ds = arceme_open(current_event, batch = batch)
+    cls = ctykeys[findfirst(ctykeys[:, 1] .== crop), 2]
+    indcrop = findall(x -> x == cls, ds.CTY.data[:, :])
+end
+
+function crop_splines(
+    current_event,
+    crop;
+    mx = 0.25,
+    nsamples = 100,
+    batch = "ARCEME-DC-8",
+    idx = "kNDVI",
+)
+    ds = arceme_open(current_event; batch)
+    event_date = arceme_eventdate(current_event)
+    r_s2 = getsplines(current_event; batch)
+
+    indcrop = getindcrop(current_event, crop; batch)
+
+    f = Figure()
+    tempo = datetime2unix.(lookup(r_s2.time))
+    tick_positions = datetime2unix.((event_date-Year(1)):Month(6):(event_date+Year(1)))  # Convert dates to Unix time (positions)
+    tick_labels = [
+        "$(Dates.format(t, "yyyy-mm-dd") )" for
+        t = (event_date-Year(1)):Month(6):(event_date+Year(1))
+    ]
+    a1 = Axis(
+        f[1, 1],
+        title = "$(current_event.uid) - $(crop)",
+        ylabel = idx,
+        xticks = (tick_positions, tick_labels),
+    )
+    seed!(42)
+    for i in rand(indcrop, nsamples)
+        ts = ds[idx].data[i[1], i[2], :]
+        mn = minimum(ts[.!isnan.(ts)])
+        lines!(
+            a1,
+            tempo,
+            r_s2[band_s2=At(idx)].data[:, i[1], i[2], 1],
+            color = mn,
+            colorrange = (0, mx),
+            colormap = (:viridis, mx),
+        )
+    end
+    Colorbar(
+        f[1, 2],
+        colorrange = (0, mx),
+        colormap = (:viridis, mx),
+        label = "Minimum $idx",
+    )
+    f
+end
+
+function tsstatsA(xout, ts, mce, mch, tempo)
+    # TO DO interpolate daily to have comparable values.
+    if (ismissing(mce) || ismissing(mch))
+        xout .= missing
+        return nothing
+    end
+    (missmask, ts) = interpolate_ts(ts, tempo)
+    ind = tempo[missmask] .> mce .&& tempo[missmask] .< mch
+    if !any(ind)
+        xout .= missing
+        return nothing
+    end
+    mx = try
+        maximum(ts[missmask][ind])
+    catch
+        @show mce
+        @show mch
+        @show tempo[missmask]
+    end
+    # tmx = Day(tempo[missmask][ind][ts[missmask][ind].==mx][1] - tempo[1]).value
+    tmx = dayofyear(tempo[missmask][ind][ts[missmask][ind] .== mx][1])
+    # need to interpolate over tempo to have a comparable basis
+
+    sts = sum(ts[missmask][ind])
+    xout .= cat(mx, tmx, sts, dims = 1)
+    return nothing
+end
+
+get_mce(ds) = map(ds.CPMCE) do x
+    (x == 0 || x >= 65500) && return missing
+    return CubePlots.yydoy2date(x)
+end
+
+get_mch(ds) = map(ds.CPMCH) do x
+    (x == 0 || x >= 65500) && return missing
+    return CubePlots.yydoy2date(x)
+end
+
+get_sce(ds) = map(ds.CPSCE) do x
+    (x == 0 || x >= 65500) && return missing
+    return CubePlots.yydoy2date(x)
+end
+
+get_val(ds; layer) = map(ds[layer]) do x
+    (x == 0 || x >= 65500) && return missing
+    return x
+end
+
+
+
+"""
+    
+"""
+function arceme_acrop(
+    current_event,
+    crop;
+    batch = "ARCEME-DC-8",
+    savef = true,
+    fname = "fig_violin_anova_$(crop)_$(current_event.uid).png",
+    printres = true,
+    io = stdout,
+)
+    ds = arceme_open(current_event; batch)
+
+    indcrop = getindcrop(current_event, crop; batch)
+    indBSB = findall(x -> x > 0 && x < 65000, ds.CPBSB.data[:, :])
+    if isempty(indcrop[indcrop .∈ (indBSB,)])
+        println(io, "No BSB data for $crop ANOVA")
+        return (nothing, nothing, nothing)
+    end
+
+    gbsb_labels = ["<35", "35-60", ">60"]
+    gbsb = map(ds.CPBSB) do x
+        (x == 0 || x >= 65500) && return missing
+        x < 35 && return 1
+        (x >= 35 || x < 60) && return 2
+        x >= 60 && return 3
+    end
+
+    mce = get_mce(ds)
+
+    mch = get_mch(ds)
+
+    outax = Dim{:mcstats}(["mx", "tmx", "sum"])
+    mcstata = xmap(
+        tsstatsA,
+        ds.kNDVI ⊘ :time_sentinel_2_l2a,
+        mce,
+        mch,
+        output = XOutput(outax, outtype = Union{Missing,Float32}),
+        function_args = (lookup(ds.time_sentinel_2_l2a),),
+    )
+    data_anova = hcat(
+        permutedims(mcstata.data[:, indcrop[indcrop .∈ (indBSB,)], 1]),
+        gbsb.data[indcrop[indcrop .∈ (indBSB,)]],
+    )
+    df_anova =
+        DataFrame(
+            mx = data_anova[:, 1],
+            tmx = data_anova[:, 2],
+            sumkndvi = data_anova[:, 3],
+            gbsb = data_anova[:, 4],
+        ) |> (df -> dropmissing!(df))
+    gdf = groupby(df_anova, :gbsb)
+    means = gdf |> (gdf -> combine(gdf, :mx => mean, :tmx => mean, :sumkndvi => mean))
+    ylabs = ["max(kNDVI)", "DOY | max(kNDVI)", "∑(kNDVI)"]
+
+    f = Figure(size = (400, 800))
+    for i = 1:3
+        ax = Axis(
+            f[i, 1],
+            xlabel = "Bare Soil Before Main Crop emergence (Days)",
+            xticks = (1:nrow(means), gbsb_labels[1:nrow(means)]),
+            ylabel = ylabs[i],
+        )
+        violin!(ax, df_anova[:, 4], df_anova[:, i], scale = :count, show_median = true)
+        for j = 1:nrow(means)
+            hlines!(
+                ax,
+                means[j, i+1],
+                xmin = (j - 1) / nrow(means),
+                xmax = j / (nrow(means)),
+                color = :grey,
+            )
+        end
+        i < 3 && hidexdecorations!(ax, grid = false)
+    end
+    Legend(
+        f[4, 1],
+        [
+            LineElement(color = :black, linestyle = nothing),
+            LineElement(color = :grey, linestyle = nothing),
+        ],
+        ["Median", "Mean"],
+        orientation = :horizontal,
+    )
+    f
+    savef && save(fname, f)
+    res = Dict()
+    for istat = 1:3
+        printres && println(io, "$(istat). $(names(df_anova)[istat]): $(ylabs[istat])\n")
+        res[names(df_anova)[istat]] = OneWayANOVATest(
+            gdf[1][!, Symbol(names(df_anova)[istat])],
+            gdf[2][!, Symbol(names(df_anova)[istat])],
+        )
+        printres && println(io, res[names(df_anova)[istat]])
+    end
+    return (f, res, data_anova)
+end
+
+function tsstatsP(xout, ts, tempo)
+    event_date = tempo[length(tempo)÷2+1]
+    (missmask, ts) = interpolate_ts(ts, tempo)
+    # missmask = findall(x -> !ismissing(x) && isfinite(x), ts)
+    indWinter = month.(tempo[missmask]) .∈ ([1, 2, 11, 12],)
+    indGrowing1 = month.(tempo[missmask]) .∈ (3:10,) .&& tempo[missmask] .< event_date
+    indGrowing2 = month.(tempo[missmask]) .∈ (3:10,) .&& tempo[missmask] .>= event_date
+    avr = mean(ts[missmask][indWinter])
+    mx1 = maximum(ts[missmask][indGrowing1])
+    tmx1 = dayofyear(tempo[missmask][indGrowing1][ts[missmask][indGrowing1] .== mx1][1])
+    sts1 = sum(ts[missmask][indGrowing1])
+    mx2 = maximum(ts[missmask][indGrowing2])
+    tmx2 = dayofyear(tempo[missmask][indGrowing2][ts[missmask][indGrowing2] .== mx2][1])
+    sts2 = sum(ts[missmask][indGrowing2])
+    xout .= cat(avr, mx1, tmx1, sts1, mx2, tmx2, sts2, dims = 1)
+    return nothing
+end
+
+
+function groupwinterkndvi(x, med)
+    x < med && return 1
+    x >= med && return 2
+end
+
+function arceme_pcrop(
+    current_event,
+    crop;
+    batch = "ARCEME-DC-6",
+    savef = true,
+    fname = "fig_violin_anova_$(crop)_$(current_event.uid).png",
+    printres = true,
+    io = stdout,
+)
+    ds = arceme_open(current_event; batch)
+    outax = Dim{:mcstats}(["avr", "mx1", "tmx1", "sum1", "mx2", "tmx2", "sum2"])
+    grapestata = xmap(
+        tsstatsP,
+        ds.kNDVI ⊘ :time_sentinel_2_l2a,
+        output = XOutput(outax, outtype = Union{Missing,Float32}),
+        function_args = (lookup(ds.time_sentinel_2_l2a),),
+    )
+    indcrop = CubePlots.getindcrop(current_event, crop; batch)
+    data_anova = permutedims(grapestata.data[:, indcrop, 1])
+    df_anova =
+        DataFrame(
+            avr = data_anova[:, 1],
+            mx1 = data_anova[:, 2],
+            tmx1 = data_anova[:, 3],
+            sumkndvi1 = data_anova[:, 4],
+            mx2 = data_anova[:, 5],
+            tmx2 = data_anova[:, 6],
+            sumkndvi2 = data_anova[:, 7],
+        ) |> (df -> dropmissing!(df))
+    med = median(df_anova.avr) # => balanced samples
+    gwinter_labels = ["< median", ">= median"]
+    transform!(df_anova, :avr => ByRow(x -> groupwinterkndvi(x, med)) => :gwinter)
+    gdf = groupby(df_anova, :gwinter)
+    # compute mean over DJF
+    means =
+        gdf |> (
+            gdf -> combine(
+                gdf,
+                :mx1 => mean,
+                :tmx1 => mean,
+                :sumkndvi1 => mean,
+                :mx2 => mean,
+                :tmx2 => mean,
+                :sumkndvi2 => mean,
+            )
+        )
+    ylabs = [
+        "max(kNDVI) | Year 1",
+        "DOY | max(kNDVI) | Year 1",
+        "∑(kNDVI) | Year 1",
+        "max(kNDVI) | Year 2",
+        "DOY | max(kNDVI) | Year 2",
+        "∑(kNDVI) | Year 2",
+    ]
+    f = Figure(size = (600, 800))
+    for i = 1:6
+        ax = Axis(
+            f[i > 3 ? i - 3 : i, i > 3 ? 2 : 1],
+            xlabel = "Average winter kNDVI",
+            xticks = (1:nrow(means), gwinter_labels[1:nrow(means)]),
+            ylabel = ylabs[i],
+        )
+        violin!(ax, df_anova[:, 8], df_anova[:, i+1], scale = :count, show_median = true)
+        for j = 1:nrow(means)
+            hlines!(
+                ax,
+                means[j, i+1],
+                xmin = (j - 1) / nrow(means),
+                xmax = j / (nrow(means)),
+                color = :grey,
+            )
+        end
+        !(i ∈ [3, 6]) && hidexdecorations!(ax, grid = false)
+    end
+    Legend(
+        f[4, 1:2],
+        [
+            LineElement(color = :black, linestyle = nothing),
+            LineElement(color = :grey, linestyle = nothing),
+        ],
+        ["Median", "Mean"],
+        orientation = :horizontal,
+    )
+    savef && save(fname, f)
+    # ANOVA
+    res = Dict()
+    for istat in names(df_anova)[2:7]
+        printres && println(io, istat)
+        res[istat] = OneWayANOVATest(gdf[1][!, Symbol(istat)], gdf[2][!, Symbol(istat)])
+        printres && println(io, res[istat])
+    end
+    return (f, res, data_anova)
+end
+
+"""
+`interpolate_ts(values::Vector, time::Vector{DateTime}; eachday=false, filter=false)`
+
+Linearly interpolates missing values in time series `values` at each time step in `time` (`eachday=false`) or at each day (`eachday=true`).
+If `filter` is `false`, values that jump for more than 0.5 over less than 20 days are discarded. If `filter` is `true`, a moving median filter is applied before interpolating (defalut `windowsize=5`)
+
+"""
+function interpolate_ts(tmp, tempo; eachday = false, filter = false, windowsize = 5)
+    dayssincestart = datetime2julian.(tempo) .- datetime2julian(tempo[1])
+    missmask = findall(x -> !ismissing(x) && isfinite(x), tmp)
+
+    if filter
+        # apply median filter
+        tmp1 = running(median, windowsize)
+        tmp[missmask] = tmp1
+    else
+        # flag jumps
+        indflag = findall(
+            broadcast(
+                (x, y) -> x > 0.5 && y <= 20,
+                diff(tmp[missmask]),
+                diff(dayssincestart[missmask]),
+            ),
+        )
+        deleteat!(missmask, indflag .+ 1)
+    end
+    # linear interpolation 
+    nodes = (dayssincestart[missmask],)
+    itp = interpolate(nodes, tmp[missmask], Gridded(Interpolations.Linear()))
+    etp = extrapolate(itp, Flat())
+    tmp1 = etp(dayssincestart)
+    return (missmask, tmp1)
+end
+
+function filter_lininterp(yout, yin, xin, nout)
+    missmask = findall(x -> !ismissing(x) && isfinite(x), yin)
+    isempty(missmask) && return yout .= missing
+    xgood = xin[missmask]
+    ygood = yin[missmask]
+    itp = interpolate((xgood,), ygood, Gridded(Interpolations.Linear()))
+    etp = extrapolate(itp, Flat())
+    yout .= etp(range(0, nout-1))
+end
+
+"""
+`plot_ts(idx::CartesianIndex, current_event::Event; layer="kNDVI")`
+
+plot time series.
+return (f, tmp1, tmp_rm, sp)
+"""
+function plot_ts(idx, current_event; layer = "kNDVI", batch = "ARCEME-DC-8")
+    ds = arceme_open(current_event, batch = batch)
+    event_date = arceme_eventdate(current_event)
+    tempo = lookup(ds.time_sentinel_2_l2a)
+    dayssincestart = datetime2julian.(tempo) .- datetime2julian(tempo[1])
+    tmp = ds[layer].data[idx[1], idx[2], :]
+    tick_positions =
+        datetime2julian.((event_date-Year(1)):Month(6):(event_date+Year(1))) .-
+        datetime2julian(tempo[1])  # Convert dates to Unix time (positions)
+    tick_labels = [
+        "$(Dates.format(t, "yyyy-mm-dd") )" for
+        t = (event_date-Year(1)):Month(6):(event_date+Year(1))
+    ]
+    f = Figure()
+    a = Axis(f[1, 1], xticks = (tick_positions, tick_labels))
+    p = scatter!(a, dayssincestart, tmp, label = "Cloud filtered")
+    missmask = findall(x -> !ismissing(x) && isfinite(x), tmp)
+    # flag jumps
+    indflag = findall(
+        broadcast(
+            (x, y) -> x > 0.5 && y <= 20,
+            diff(tmp[missmask]),
+            diff(dayssincestart[missmask]),
+        ),
+    )
+    deleteat!(missmask, indflag .+ 1)
+    # scatterlines!(a, dayssincestart[missmask], tmp[missmask], label = "Jump filtered")
+    # linear interpolation 
+    tmp1 = fill(NaN, length(tmp))
+    tmp1[missmask] .= tmp[missmask]
+    nodes = (dayssincestart[missmask],)
+    itp = interpolate(nodes, tmp[missmask], Gridded(Interpolations.Linear()))
+    if missmask[1] > 1
+        tmp1[1:(missmask[1]-1)] .= tmp[missmask[1]]
+    end
+    if missmask[end] < length(dayssincestart)
+        tmp1[(missmask[end]+1):end] .= tmp[missmask[end]]
+    end
+    tmp1[missmask[1]:missmask[end]] = itp(dayssincestart[missmask[1]:missmask[end]])
+
+    scatterlines!(
+        a,
+        dayssincestart,
+        tmp1,
+        color = Makie.wong_colors()[1],
+        markersize = 5,
+        label = "Linear interpolation",
+    )
+    f
+    # # 2. SavitzkyGolay filter 
+    # tmp_sg1 = savitzky_golay(tmp1, 7, 2)
+    # lines!(a, dayssincestart, tmp_sg1.y, label = "Savitzky–Golay w=3, o=1")
+    # # tmp_sg2 = savitzky_golay(tmp_sg1.y, 5, 3)
+    # # lines!(dayssincestart, tmp_sg2.y, label = "Savitzky–Golay w=5, o=3")
+    # # tmp_sg3 = savitzky_golay(tmp_sg2.y, 7, 2)
+    # # lines!(dayssincestart, tmp_sg3.y, label = "Savitzky–Golay w=7, o=5")
+    # # tmp_sg4 = savitzky_golay(tmp_sg3.y, 9, 6)
+    # # lines!(dayssincestart, tmp_sg4.y, label = "Savitzky–Golay w=9, o=6")
+    # f
+    # # the succesive iterations don't do anything! Because of the linear interpolation, I suppose?
+    # # Maybe I should use a rolling mean before
+
+    # running mean
+    tmp_rm = running(mean, tmp1, 5)
+    lines!(
+        a,
+        dayssincestart[1:(end-2)],
+        tmp_rm[3:end],
+        color = Makie.wong_colors()[2],
+        label = "Rolling mean w=5",
+    )
+
+    # running median before interpolation
+    tmp_rmd = fill(NaN, length(tmp))
+    w=3
+    tmp_rmd[missmask] = running(median, tmp[missmask], w)
+
+    itp = interpolate(nodes, tmp_rmd[missmask], Gridded(Interpolations.Linear()))
+    # pad with first value
+    if missmask[1] > 1
+        tmp_rmd[1:(missmask[1]-1)] .= tmp_rmd[missmask[1]]
+    end
+    # pad with last value
+    if missmask[end] < length(dayssincestart)
+        tmp_rmd[(missmask[end]+1):end] .= tmp_rmd[missmask[end]]
+    end
+    tmp_rmd[missmask[1]:missmask[end]] = itp(dayssincestart[missmask[1]:missmask[end]])
+    scatterlines!(
+        a,
+        dayssincestart[missmask],
+        tmp_rmd[missmask],
+        color = Makie.wong_colors()[4],
+        label = "running median w=$w",
+    )
+    scatterlines!(
+        a,
+        dayssincestart,#[1:(end-2)],
+        tmp_rmd,#[3:end],
+        color = Makie.wong_colors()[4],
+        markersize = 5,
+        label = "interpolated running median w=$w",
+    )
+
+    # running median after interpolarion
+    tmp_rmd2 = running(median, tmp1, 5)
+    lines!(
+        a,
+        dayssincestart[1:(end-2)],
+        tmp_rmd2[3:end],
+        color = Makie.wong_colors()[5],
+        label = "Rolling median w=5",
+    )
+
+    r_s2 = CubePlots.getsplines(current_event)
+    sp = r_s2[band_s2=At(layer)].data[:, idx[1], idx[2], 1]
+    lines!(a, sp, color = Makie.wong_colors()[3], label = "Spline 19 knots")
+    Legend(f[2, 1], a, orientation = :horizontal, nbanks = 3)
+    f
+    return (f, tmp1, tmp_rm, sp)
+end
+
+export plotlc,
+    plotcty,
+    plotmceh,
+    piecty,
+    fpfig,
+    getsplines,
+    plot_cpbsb_splines,
+    crop_splines,
+    getindcrop,
+    arceme_acrop,
+    arceme_pcrop,
+    plot_ts
+end
\ No newline at end of file
diff --git a/src/plots/fingerprint_plots.jl b/src/plots/fingerprint_plots.jl
index 6202cdf..18a1369 100644
--- a/src/plots/fingerprint_plots.jl
+++ b/src/plots/fingerprint_plots.jl
@@ -9,6 +9,7 @@ using ..ARCEMEAnalysis: local_cubepath, arceme_legends, Event
 using DataStructures: counter, OrderedDict
 using Statistics: mean
 using GeoMakie
+using DiskArrays: cache
 import DimensionalData as DD
 
 function filter_spline(xout, samplets, x, nout)
@@ -33,94 +34,92 @@ end
 
 
 
-
-pairnow = arceme_validpairs()[4]
-
-
-
-function most_common_lc(pairnow)
+function most_common_class(pairnow)
     ds1, ds2 = arceme_open.(pairnow)
     #Find maximum landcover in both cubes
-    lcfrac1, lcfrac2 = (ds1.lc_fraction[:], ds2.lc_fraction[:])
-    (_, i1), (_, i2) = findmax.((lcfrac1, lcfrac2))
-    if lcfrac1[i1] * lcfrac2[i1] > lcfrac1[i2] * lcfrac2[i2]
-        ds1.lc.val[i1]
+    classfrac1, classfrac2 = (ds1.class_fractions[:], ds2.class_fractions[:])
+    (_, i1), (_, i2) = findmax.((classfrac1, classfrac2))
+    if classfrac1[i1] * classfrac2[i1] > classfrac1[i2] * classfrac2[i2]
+        ds1.class.val[i1]
     else
-        ds1.lc.val[i2]
+        ds1.class.val[i2]
     end
 end
 
 function pairwise_without_missings(x, y)
     out = Point2f[]
+    indices = Int[]
     for i in eachindex(x, y)
         if !ismissing(y[i]) && isfinite(y[i])
             push!(out, Point2f(x[i], y[i]))
+            push!(indices, i)
         end
     end
-    out
+    out, indices
 end
 
-function open_plot_data(ev, lc; strata="ESA_LC")
-    @info "Opening $ev"
-    ds = arceme_open(ev)
-    firstday = arceme_eventdate(ev) - Year(1)
-    lastday = arceme_eventdate(ev) + Year(1)
+function open_plot_data(ds, class; strata="ESA_LC")
+    firstday = arceme_eventdate(ds) - Year(1)
+    lastday = arceme_eventdate(ds) + Year(1)
     nms = Millisecond(lastday - firstday).value
     linscal = d -> Millisecond(d - firstday).value * 2 / nms
     x_s2 = map(linscal, ds.time_sentinel_2_l2a.val)
     x_s1 = map(linscal, ds.time_sentinel_1_rtc.val)
     outax = YAXArrays.time(range(firstday, lastday, step=Day(1)))
 
-    s2 = ds.s2_indices[lc=DD.At(lc)]
-    s1 = ds.s1_indices[lc=DD.At(lc)]
+    s2 = ds.s2_indices[class=DD.At(class)]
+    s1 = ds.s1_indices[class=DD.At(class)]
 
     r_s2 = xmap(filter_spline, s2 ⊘ :time_sentinel_2_l2a, output=XOutput(outax), function_args=(x_s2, length(outax),))
     r_s1 = xmap(filter_spline, s1 ⊘ :time_sentinel_1_rtc, output=XOutput(outax), function_args=(x_s1, length(outax),))
 
     dataplots = OrderedDict()
     interplots = OrderedDict()
+    index_orig = OrderedDict()
 
     for band in ds.band_s2
-
-        interplots[band] = pairwise_without_missings(range(0, 2, length=length(outax)), r_s2[band_s2=DD.At(band)].data[:])
-        dataplots[band] = pairwise_without_missings(x_s2, s2[band_s2=DD.At(band)].data[:])
+        interplots[band], _ = pairwise_without_missings(range(0, 2, length=length(outax)), r_s2[band_s2=DD.At(band)].data[:])
+        dataplots[band], index_orig[band] = pairwise_without_missings(x_s2, s2[band_s2=DD.At(band)].data[:])
     end
 
     for band in ds.band_s1
-        interplots[band] = pairwise_without_missings(range(0, 2, length=length(outax)), r_s1[band_s1=DD.At(band)].data[:])
-        dataplots[band] = pairwise_without_missings(x_s1, s1[band_s1=DD.At(band)].data[:])
+        interplots[band], _ = pairwise_without_missings(range(0, 2, length=length(outax)), r_s1[band_s1=DD.At(band)].data[:])
+        dataplots[band], index_orig[band] = pairwise_without_missings(x_s1, s1[band_s1=DD.At(band)].data[:])
     end
 
-    lcfracs = ds.lc_fraction.data[:]
-    sfracs = sortperm(lcfracs, rev=true)
+    classfracs = ds.class_fractions.data[:]
+    sfracs = sortperm(classfracs, rev=true)
     strs = map(sfracs[1:3]) do ifrac
-        string(collect(values(arceme_legends[strata]))[ifrac], " => ", round(lcfracs[ifrac] * 100), "%")
+        string(collect(values(arceme_legends[strata]))[ifrac], " => ", round(classfracs[ifrac] * 100), "%")
     end
-    lcdesc = join(strs, ", ")
+    classdesc = join(strs, ", ")
 
 
-    (dataplots, interplots, lcdesc)
+    (dataplots, interplots, classdesc, index_orig)
 end
 
 
-function fingerprint_plot(pairid=1; lc=nothing, plotdict=IdDict(), interactive=true, strata="ESA_LC")
+function fingerprint_plot(pairid=1; class=nothing, plotdict=IdDict(), interactive=true, strata="ESA_LC", mask_clouds=false)
 
     ipair = Observable{Int}(pairid)
     alleventpairs = arceme_validpairs()
 
-    if lc === nothing
-        lc = most_common_lc(alleventpairs[pairid])
+    if class === nothing
+        class = most_common_class(alleventpairs[pairid])
     end
 
-    lc_choos = Observable{Any}(lc)
+    class_choos = Observable{Any}(class)
+
+    ds = lift(ipair) do i
+        ev = alleventpairs[i]
+        cache.(arceme_open.(ev))
+    end
 
-    d = lift(ipair, lc_choos) do i, lc
-        ev = alleventpairs[i][1]
-        open_plot_data(ev, lc)
+    d = lift(ds, class_choos) do dss, class
+        open_plot_data(dss[1], class)
     end
-    dhp = lift(ipair, lc_choos) do i, lc
-        ev = alleventpairs[i][2]
-        open_plot_data(ev, lc)
+    dhp = lift(ds, class_choos) do dss, class
+        open_plot_data(dss[2], class)
     end
     fig = Figure(size=(1600, 2600))
 
@@ -145,15 +144,19 @@ function fingerprint_plot(pairid=1; lc=nothing, plotdict=IdDict(), interactive=t
         end
     end
 
+    repr_ind_d = Observable{Any}(nothing)
+    repr_ind_dhp = Observable{Any}(nothing)
+    omask_clouds = Observable{Bool}(mask_clouds)
+
     npl = maximum(first, keys(allax))
     for ipl in 1:npl
         linkyaxes!(allax[(ipl, 1)], allax[(ipl, 2)])
     end
     ov1 = Axis(fig[5:6, 3], aspect=DataAspect(), title="Drought")
-    img1 = lift((y, i) -> arceme_representative_image(alleventpairs[i][1], lc=y), lc_choos, ipair)
+    img1 = lift((y, i, ind, mask_clouds) -> arceme_representative_image(i[1]; class=y, repr_ind=ind, mask_clouds), class_choos, ds, repr_ind_d, omask_clouds)
     heatmap!(ov1, img1)
     ov2 = Axis(fig[7:8, 3], aspect=DataAspect(), title="Drought + HP")
-    img2 = lift((y, i) -> arceme_representative_image(alleventpairs[i][2], lc=y), lc_choos, ipair)
+    img2 = lift((y, i, ind, mask_clouds) -> arceme_representative_image(i[2]; class=y, repr_ind=ind, mask_clouds), class_choos, ds, repr_ind_dhp, omask_clouds)
     heatmap!(ov2, img2)
 
 
@@ -167,18 +170,38 @@ function fingerprint_plot(pairid=1; lc=nothing, plotdict=IdDict(), interactive=t
     # xlims!(gax, lift(xpos -> minimum(xpos) - 20, xpos), lift(xpos -> maximum(xpos) + 20, xpos))
     # ylims!(gax, lift(ypos -> minimum(ypos) - 20, ypos), lift(ypos -> maximum(ypos) + 20, ypos))
     if interactive
-        menupair = Textbox(fig, placeholder=string(ipair[]), validator=Int, tellwidth=false)
+
         #menupair = Slider(fig, range=1:length(alleventpairs), startvalue=pairid, update_while_dragging=false)
-        menulc = Menu(fig, options=collect(values(arceme_legends[strata])), default=lc)
+
+        menupair = Textbox(fig, placeholder=string(ipair[]), validator=Int)
+        menuclass = Menu(fig, options=collect(values(arceme_legends[strata])), default=class)
+        # cb = Checkbox(checkboxes[1, 1], checked=false)
+        cb2 = Checkbox(fig, checked=true)
+        # cb3 = Checkbox(checkboxes[3, 1], checked=true)
+        on(c->setindex!(omask_clouds, c), cb2.checked)
+
         fig[1:2, 3] = vgrid!(
             menupair,
             Label(fig, lift(i -> i[3], d), color=:red),
             Label(fig, lift(i -> i[3], dhp), color=:blue),
-            menulc,
-
+            menuclass,
+            cb2,
         )
-        on(menulc.selection) do s
-            lc_choos[] = s
+
+
+        #fig[1:2,3] = allmenus
+
+
+        # checkboxes = GridLayout(allmenus[5:7, 1:3])
+
+
+
+
+        # Label(checkboxes[1, 2], "Classmask", halign=:left)
+        # Label(checkboxes[2, 2], "Cloudmask", halign=:left)
+        # Label(checkboxes[3, 2], "SCL Mask", halign=:left)
+        on(menuclass.selection) do s
+            class_choos[] = s
         end
         on(menupair.stored_string) do v
             ipair[] = parse(Int, v)
@@ -188,47 +211,77 @@ function fingerprint_plot(pairid=1; lc=nothing, plotdict=IdDict(), interactive=t
             autolimits!(ov1)
             autolimits!(ov2)
         end
+        on(events(fig).mousebutton) do event
+            if event.button == Mouse.left && event.action == Mouse.press
+                p, i = pick(fig)
+                index, i_dhp = get(plotdict, p, (nothing, nothing))
+                if !isnothing(index)
+                    if i_dhp == 1
+                        indices = last(d[])
+                        repr_ind_d[] = indices[index][i]
+                    else
+                        indices = last(dhp[])
+                        repr_ind_dhp[] = indices[index][i]
+                    end
+                    Consume(true)
+                end
+            end
+        end
     end
     fig
 end
 
 
-arceme_representative_image(ev::Event; index_use="NDVI", brighten_factor=2, lc=nothing) =
-    arceme_representative_image(arceme_open(ev); index_use, brighten_factor, lc)
+arceme_representative_image(ev::Event; index_use="NDVI", brighten_factor=2, repr_ind=nothing, class=nothing, mask_clouds=false) =
+    arceme_representative_image(arceme_open(ev); index_use, brighten_factor, repr_ind, class, mask_clouds)
 
 """
     arceme_representative_image(ev)
 
 Creates an RGB image highlighting on a cloud-free time step with high NDVI for a certain land cover type 
 """
-function arceme_representative_image(ds; index_use="NDVI", brighten_factor=2, lc=nothing, strata="ESA_LC")
-    if lc === nothing
-        lcfrac = ds.lc_fraction.data[:]
-        _, ilc = findmax(lcfrac)
-        lc = ds.lc.val[ilc]
+function arceme_representative_image(ds; index_use="NDVI", brighten_factor=2, repr_ind=nothing, class=nothing, strata="ESA_LC", mask_clouds=false)
+    if class === nothing
+        classfrac = ds.class_fractions.data[:]
+        _, iclass = findmax(classfrac)
+        class = ds.class.val[iclass]
     end
-    ndvi = ds.s2_indices[band=DD.At(index_use), lc=DD.At(lc)]
-    allndviranks = sortperm(ndvi.data[:], rev=true)
-    repr_ind = 0
-    cloudfrac = ds.cloud_fraction.data[:]
-    for i in 1:length(allndviranks)
-        icandidate = findfirst(==(i), allndviranks)
-        if !ismissing(ndvi[icandidate]) && cloudfrac[icandidate] < 0.05
-            repr_ind = icandidate
-            break
+    if repr_ind === nothing
+        ndvi = ds.s2_indices[band=DD.At(index_use), class=DD.At(class)]
+        allndviranks = sortperm(ndvi.data[:], rev=true)
+        repr_ind = 0
+        cloudfrac = ds.cloud_fraction.data[:]
+        for i in 1:length(allndviranks)
+            icandidate = findfirst(==(i), allndviranks)
+            if !ismissing(ndvi[icandidate]) && cloudfrac[icandidate] < 0.05
+                repr_ind = icandidate
+                break
+            end
         end
-    end
-    if repr_ind == 0
-        return YAXArray((ds.x, ds.y), fill(RGBA(0.0, 0.0, 0.0, 1.0), 1000, 1000), Dict{Any,Any}())
+        if repr_ind == 0
+            return YAXArray((ds.x, ds.y), fill(RGBA(0.0, 0.0, 0.0, 1.0), 1000, 1000), Dict{Any,Any}())
+        end
+    else
+        @show repr_ind
+        @show ds.time_sentinel_2_l2a[repr_ind]
     end
     brfilt(a, factor, alpha) = RGBA(clamp(a.r * factor, 0, 1), clamp(a.g * factor, 0, 1), clamp(a.b * factor, 0, 1), alpha)
-    ilc = findfirst(==(lc), collect(values(arceme_legends[strata])))
-    r = broadcast(arceme_rgb(ds)[time_sentinel_2_l2a=repr_ind], ds.ESA_LC[:, :, 1], brighten_factor) do col, lc, f
+    ilc = findfirst(==(class), collect(values(arceme_legends[strata])))
+    r = broadcast(arceme_rgb(ds)[time_sentinel_2_l2a=repr_ind], ds[strata][:, :, 1], brighten_factor) do col, lc, f
         lcv = ARCEMEAnalysis.lckeymap(lc)
         alpha = lcv == ilc ? 1.0 : 0.5
         brfilt(col, f, alpha)
     end
-    r[:, :]
+    if mask_clouds
+        cl = ds.cloud_mask[:, :, repr_ind]
+        scl = ds.SCL[:, :, repr_ind]
+        is_cloud = map(ARCEMEAnalysis._is_cloud, cl.data, scl.data)
+        res = r[:, :]
+        res[is_cloud] .= RGBA(0.0, 0.0, 1.0, 1.0)
+        res
+    else
+        r[:, :]
+    end
 end
 
 function compute_knots(xgood, ntarget=19)
diff --git a/src/s1_helpers.jl b/src/s1_helpers.jl
index b453071..281162c 100644
--- a/src/s1_helpers.jl
+++ b/src/s1_helpers.jl
@@ -7,33 +7,34 @@ Analyses the sentinel 1 time stamps of the datasets and groups them into series
 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)
+function arceme_s1_position(ds; distance_threshold=61)
     ts = ds.time_sentinel_1_rtc.val
-    sixdays = 6*24*60*60*1000
+    sixdays = 6 * 24 * 60 * 60 * 1000
     group_offsets = [0]
     groups = [[1]]
     for its in 2:length(ts)
         groupfound = false
+        fac = NaN
         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
+            fac = (Millisecond(ts[its] - ts[1]).value - group_offsets[igroup]) / sixdays
+            if abs(round(fac) - fac) < (distance_threshold * 1000 / sixdays)
+                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])
+            newoffset = mod(Millisecond(ts[its] - ts[1]).value, sixdays) #- floor(Int, fac) * sixdays
+            push!(group_offsets, newoffset)
+            push!(groups, [its])
         end
     end
-    positions = zeros(Int,length(ts))
+    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)
+    ds.cubes[:s1_position] = YAXArray((ds.time_sentinel_1_rtc,), positions)
     nothing
 end
 
diff --git a/src/spatialdebias.jl b/src/spatialdebias.jl
index c08c1d6..f8e5aa3 100644
--- a/src/spatialdebias.jl
+++ b/src/spatialdebias.jl
@@ -1,28 +1,29 @@
 export arceme_bias_corrected_fp, arceme_uncorrected_fp
 
-function sincosparams(c,timname = :time_sentinel_2_l2a)
-    fitmat = bare_matrix(c,timname)
+function sincosparams(c, timname=:time_sentinel_2_l2a)
+    fitmat = bare_matrix(c, timname)
     xmap(
-        fit_sincos, 
-        c ⊘ timname, 
-        output=XOutput(DD.Dim{:param}(["p1","p2","p3"]),outtype=Float32),
+        fit_sincos,
+        c ⊘ timname,
+        output=XOutput(DD.Dim{:param}(["p1", "p2", "p3"]), outtype=Float32),
         inplace=false,
-        function_args=(fitmat,)
+        function_args=(fitmat,),
+        allow_threads=true
     )
 end
 
-function fit_sincos(ts,fitmat)
-    igood = findall(!isnan,ts)
-    view(fitmat,igood,:) \ view(ts,igood)
+function fit_sincos(ts, fitmat)
+    igood = findall(!isnan, ts)
+    view(fitmat, igood, :) \ view(ts, igood)
 end
 
-function bare_matrix(c,timname)
+function bare_matrix(c, timname)
     ms_per_year = 365.25*24*60*60*1000
-    timdim = DD.dims(c,timname)
+    timdim = DD.dims(c, timname)
     t = map(timdim.val) do x
         (x-first(timdim)).value/ms_per_year*2pi
     end
-    [sin.(t) cos.(t) fill(1.0,length(t))]
+    [sin.(t) cos.(t) fill(1.0, length(t))]
 end
 
 function nstrata(strata)
@@ -40,51 +41,64 @@ end
 
 Computes a cloud-biased corrected footprint of `dataset[band]` aggregated by stratification band `strata`. 
 """
-function arceme_bias_corrected_fp(band::String, dataset::Dataset; strata="ESA_LC", timeaxis=:time_sentinel_2_l2a)
+function arceme_bias_corrected_fp(band::String, dataset::Dataset; strata="ESA_LC", timeaxis=:time_sentinel_2_l2a, fpstat=mean)
     ncl = nstrata(strata)
     lccube = lckeymap(dataset, strata=strata)
-    cloudcube = dataset.cloud_mask
-    sclcube = dataset.SCL
+    cloudhalocube = dataset.cloudmask_halo
     inputcube = dataset[band]
-    pars = ARCEMEAnalysis.sincosparams(inputcube)[:,:,:,1]
-    timdim = DD.dims(inputcube,timeaxis)
-    fitmat = YAXArray((timdim, pars.param),ARCEMEAnalysis.bare_matrix(inputcube,timeaxis))
-
-    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=1:ncl)
-    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 _is_cloud(cl,scl)
-            return NaN
+    pars = ARCEMEAnalysis.sincosparams(inputcube)[:, :, :, 1]
+    timdim = DD.dims(inputcube, timeaxis)
+    fitmat = YAXArray((timdim, pars.param), ARCEMEAnalysis.bare_matrix(inputcube, timeaxis))
+
+    # 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=1:ncl)
+    # 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 _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=1:ncl)
+    # fp_clouded_expected = YAXArrays.compute(mean.(win_clouded).data)[:, 1, 1, 1, :]
+    varax = DD.Dim{:Variable}(["expected", "expected_cloudfiltered", "data_cloudfiltered"])
+    predcube = xmap(pars ⊘ :param, fitmat ⊘ :param, inputcube, cloudhalocube, inplace=true, output=XOutput(varax, outtype=Float32)) do out, p, t, x, clh
+        expected = p[1] * t[1] + p[2] * t[2] + p[3] * t[3]
+        out[1] = expected
+        if clh != 1
+            out[2] = NaN
+            out[3] = NaN
         else
-            p[1]*t[1]+p[2]*t[2]+p[3]*t[3]
+            out[2] = expected
+            out[3] = x
         end
     end
-    win_clouded = YAXArrays.windows(clouded_expected, lccube, expected_groups=1:ncl)
-    fp_clouded_expected = mean.(win_clouded)[:,1,1,1,:].data
+    win_data = YAXArrays.windows(predcube, lccube, expected_groups=1:ncl)
+    fp = YAXArrays.compute(xmap(mean, win_data, min_nvalid=50000).data, showprogress=false)[:, :, 1, 1, 1, :]
+    # win_data = YAXArrays.windows(inputcube_filtered, lccube, expected_groups=1:ncl)
+    # fp = YAXArrays.compute(mean.(win_data).data)[:, 1, 1, :]
 
-    inputcube_filtered = xmap(inputcube, cloudcube,sclcube,inplace=false,output=XOutput(outtype=Float32)) do x,cl,scl
-        _is_cloud(cl,scl) ? NaN : x
-    end
-    win_data = YAXArrays.windows(inputcube_filtered, lccube, expected_groups=1:ncl)
-    fp = mean.(win_data)[:,1,1,:].data
 
-    newdata = fp[:,:] .+ fp_clearsky_expected[:,:] .- fp_clouded_expected[:,:]
+    newdata = fp[:, 3, :] .+ fp[:, 1, :] .- fp[:, 2, :]
 
     classax = classaxis(strata, ncl)
 
+    abundance = counter(lccube)
+
     Dataset(
-        fp = YAXArray((classax,timdim),newdata),
-        fp_uncorrected=YAXArray((classax, timdim), fp[:, :]),
-        fp_clearsky_expected=YAXArray((classax, timdim), fp_clearsky_expected[:, :]),
-        fp_clouded_expected=YAXArray((classax, timdim), fp_clouded_expected[:, :]),
+        fp=YAXArray((classax, timdim), newdata),
+        fp_uncorrected=YAXArray((classax, timdim), fp[:, 3, :]),
+        fp_clearsky_expected=YAXArray((classax, timdim), fp[:, 1, :]),
+        fp_clouded_expected=YAXArray((classax, timdim), fp[:, 2, :]),
         params=pars,
         smooth_matrix=fitmat,
+        class_fractions=YAXArray((classax,), [abundance[i] for i in 1:ncl])
     )
-    
+
 end
 
 """
@@ -92,14 +106,14 @@ end
 
 Computes the footprint of `dataset[band]` aggregated by stratification classes in band `strata`. 
 """
-function arceme_uncorrected_fp(band, dataset; strata="ESA_LC", timeaxis=:time_sentinel_1_rtc)
+function arceme_uncorrected_fp(band, dataset; strata="ESA_LC", timeaxis=:time_sentinel_1_rtc, fpstat=mean)
     ncl = nstrata(strata)
     lccube = lckeymap(dataset, strata=strata)
     inputcube = dataset[band]
-    timdim = DD.dims(inputcube,timeaxis)
-    
+    timdim = DD.dims(inputcube, timeaxis)
+
     win = YAXArrays.windows(inputcube, lccube, expected_groups=1:ncl)
-    fp = mean.(win)[:,1,1,:].data
+    fp = fpstat.(win)[:, 1, 1, :].data
 
     abundance = counter(lccube)
     classax = classaxis(strata, ncl)
@@ -109,5 +123,5 @@ function arceme_uncorrected_fp(band, dataset; strata="ESA_LC", timeaxis=:time_se
         fp=YAXArray((classax, timdim), fp),
         class_fractions=YAXArray((classax,), [abundance[i] for i in 1:ncl])
     )
-    
+
 end
\ No newline at end of file
diff --git a/src/spectral_helpers.jl b/src/spectral_helpers.jl
index ba2fc01..395f93f 100644
--- a/src/spectral_helpers.jl
+++ b/src/spectral_helpers.jl
@@ -16,17 +16,19 @@ function (ntw::NTWrapper{F,names})(args...) where {F,names}
 end
 NTWrapper(f,names::Tuple,consts::NamedTuple,indices) = NTWrapper(f,Val(names),consts,indices)
 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...)
+    _is_clear_or_shadow(bands.cl, bands.scl) || return map(_ -> NaN, indices_tuple)
+    bandparams = map(ARCEMEAnalysis.boa, bands)
+    bands2 = map(between_one_and_zero, bandparams)
+    allparams = (; bands2..., consts...)
     listofindices(indices_tuple, allparams)
 end
+between_one_and_zero(x) = ifelse(0.0 <= x <= 1.0, x, oftype(x, NaN))
 function inner_compute_indices_s1(bands, consts, indices_tuple)
     allparams = (; bands..., consts...)
     listofindices(indices_tuple, allparams)
 end
 
-function _compute_indices(ds,indices,platform)
+function _compute_indices(ds, indices, platform; showprog=true)
     indices_tuple = ((SI.indices[k] for k in indices)...,)
 
     allvars = mapreduce(r->Set(string.(SI._band_names(r))),union!,indices_tuple)
@@ -50,7 +52,9 @@ function _compute_indices(ds,indices,platform)
 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 compute_indexx(index::SpectralIndices.SpectralIndex{<:Any,B}, params::NamedTuple) where B
+    index.compute(Float64, params[B]...)
+end
 function listofindices(indices, values)
     map(indices) do index
         compute_indexx(index, values)
@@ -76,7 +80,7 @@ function estimate_sigma_per_lc(ds;bands=(:B08,:B04))
     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 = DAE.compute(DAE.results_as_diskarrays(op)[1], showprogress=false)
     res[1, 1, 1, :]
 end