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/Project.toml b/Project.toml
index f0bfccd..fff9fe6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,7 +8,9 @@ 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"
@@ -18,6 +20,8 @@ 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"
@@ -26,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"
@@ -39,7 +45,9 @@ 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"
@@ -49,6 +57,8 @@ 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"
@@ -57,7 +67,9 @@ 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"
diff --git a/src/ARCEMEAnalysis.jl b/src/ARCEMEAnalysis.jl
index f8b735f..b45e995 100644
--- a/src/ARCEMEAnalysis.jl
+++ b/src/ARCEMEAnalysis.jl
@@ -1,12 +1,22 @@
 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
@@ -17,8 +27,14 @@ 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)
@@ -26,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",
@@ -53,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(
@@ -98,26 +107,26 @@ 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) = ifelse(k > 90, (k + 10), k) ÷ 10 + 1
-function lckeymap(ds::Dataset; strata="ESA_LC")
+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)
@@ -127,19 +136,38 @@ function lckeymap(ds::Dataset; strata="ESA_LC")
 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, arceme_cloudmask_halo, arceme_precompute_fingerprints
+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-DHP-GLOBAL")
+function _arceme_cubenames(; batch = "ARCEME-DC-DHP-GLOBAL")
     cubenames = if local_cubepath === nothing
         store = S3Store("$(batch)/", MinioConfig(httpstore))
 
@@ -183,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,
@@ -203,7 +233,10 @@ function arceme_eventlist(; csv=joinpath(@__DIR__, "..", "data", "dhp_global_sub
     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,
@@ -211,7 +244,8 @@ function arceme_eventlist(source::String; geojson=joinpath(@__DIR__, "..", "data
             row.geometry[1],
             row.geometry[2],
             DateTime(row.startdate),
-            Symbol(source),)
+            Symbol(source),
+        )
     end
 end
 
@@ -224,7 +258,8 @@ 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-DHP-GLOBAL")
@@ -232,25 +267,56 @@ arceme_eventpairs() =
 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-DHP-GLOBAL")
+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-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)
+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-DHP-GLOBAL")
+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-DHP-GLOBAL") = open_dataset(joinpath(local_cubepath, "$batch-fingerprints.zip"), path=cubename, force_datetime=true)
+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)
@@ -261,10 +327,16 @@ 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-DHP-GLOBAL") = arceme_landcover(arceme_open(ev, batch=batch))
+arceme_landcover(ev::Event; batch = "ARCEME-DC-DHP-GLOBAL") =
+    arceme_landcover(arceme_open(ev, batch = batch))
 
 
 """ 
@@ -277,11 +349,21 @@ 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-DHP-GLOBAL", 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")), force_datetime=true)
+        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
@@ -289,15 +371,37 @@ function arceme_open(cubename; batch="ARCEME-DC-DHP-GLOBAL", indices=true, trylo
             #     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)
+                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]
+                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)
+        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
@@ -305,13 +409,19 @@ function arceme_open(cubename; batch="ARCEME-DC-DHP-GLOBAL", indices=true, trylo
         if fingerprint
             try
                 ds_fp = arceme_open_fingerprint(cubename; batch)
-                main_ds = Dataset(; main_ds.cubes..., ds_fp.cubes..., properties=main_ds.properties)
+                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
@@ -354,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
 
@@ -380,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)
@@ -406,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)
@@ -438,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
 
 """
@@ -476,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
 
 
@@ -491,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...)
@@ -514,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
 
 """
@@ -535,13 +694,18 @@ 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-DHP-GLOBAL", subset=:)
-    @showprogress pmap(arceme_validpairs(batch=batch)[subset]) do ev
-        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
             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_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)
@@ -564,20 +728,28 @@ 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, showprogress=false, use_dagger=false)
+    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)
+        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")))
+    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)
+    rm(joinpath(output_base, name), recursive = true)
 end
 
 # Utility function to zip a directory
-function apply_dirfiles(f, basepath, dir="")
+function apply_dirfiles(f, basepath, dir = "")
     flist = readdir(joinpath(basepath, dir))
     for file in flist
         if isdir(joinpath(basepath, dir, file))
@@ -601,7 +773,12 @@ end
 
 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-DHP-GLOBAL")
+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
@@ -609,12 +786,13 @@ function arceme_create_indexcubes(event_list; indices_s1=["DpRVIVV"], indices_s2
 
     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)
@@ -627,11 +805,17 @@ 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
@@ -640,21 +824,37 @@ function arceme_index_fingerprints(ds; indices_s1=["DpRVIVV", "vv_db", "vh_db"],
     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)
+            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)
+        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
     if !isempty(indices_s2)
         indexdim_s2 = DD.Dim{:band_s2}(indices_s2)
 
         fp_s2 = map(indices_s2) do band
-            arceme_bias_corrected_fp(band, ds, strata=strata)
+            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)))
+        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)
+        s2cubes = (
+            :s2_indices => fingerprint_sparse_s2,
+            :uncorrected_s2_indices => fingerprint_uncor_sparse_s2,
+            :class_fractions => fp_s2[1].class_fractions,
+        )
     end
     Dataset(; s1cubes..., s2cubes...)
 end
@@ -663,7 +863,8 @@ 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"]))
+    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
@@ -672,7 +873,12 @@ function arceme_fractions(ds)
     end
     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
+    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)
@@ -684,22 +890,52 @@ 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
+    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)),)
+    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"))
+    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")
+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)
+        savedataset(fp, path = joinpath(tmpdir, arceme_cubename(ev)), overwrite = true)
     end
 
     outputzipfile = joinpath(ARCEMEAnalysis.local_cubepath, "$batch-fingerprints.zip")
@@ -713,9 +949,11 @@ 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/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/spatialdebias.jl b/src/spatialdebias.jl
index ff19ce1..f8e5aa3 100644
--- a/src/spatialdebias.jl
+++ b/src/spatialdebias.jl
@@ -1,29 +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,),
         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)
@@ -41,20 +41,20 @@ 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)
     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))
+    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, lazy=true) do p, t
+    # 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 = YAXArrays.compute(mean.(win_clearsky).data)[:, 1, 1, 1, :]
+    # 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)
@@ -77,12 +77,12 @@ function arceme_bias_corrected_fp(band::String, dataset::Dataset; strata="ESA_LC
             out[3] = x
         end
     end
-
     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, :]
 
+
     newdata = fp[:, 3, :] .+ fp[:, 1, :] .- fp[:, 2, :]
 
     classax = classaxis(strata, ncl)
@@ -90,7 +90,7 @@ function arceme_bias_corrected_fp(band::String, dataset::Dataset; strata="ESA_LC
     abundance = counter(lccube)
 
     Dataset(
-        fp = YAXArray((classax,timdim),newdata),
+        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, :]),
@@ -98,7 +98,7 @@ function arceme_bias_corrected_fp(band::String, dataset::Dataset; strata="ESA_LC
         smooth_matrix=fitmat,
         class_fractions=YAXArray((classax,), [abundance[i] for i in 1:ncl])
     )
-    
+
 end
 
 """
@@ -106,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)
@@ -123,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