Melwey

Project.toml

@@ -9,12 +9,15 @@ Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
+DiskArrayEngine = "2d4b2e14-ccd6-4284-b8b0-2378ace7c126"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Minio = "4281f0d9-7ae0-406e-9172-b7277c1efa20"
+OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
+SpectralIndices = "df0093a1-273d-40bc-819a-796ec3476907"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StringViews = "354b36f9-a18e-4713-926e-db85100087ba"
 YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c"
@@ -24,18 +27,21 @@ ZipArchives = "49080126-0e18-4c2a-b176-c102e4b3760c"
 [compat]
 CSV = "0.10.15"
 Colors = "0.13.1"
-DataStructures = "0.18.22"
+DataStructures = "0.18.22, 0.19"
 Dates = "1.11.0"
 DimensionalData = "0.29.24"
+DiskArrayEngine = "0.3.0"
 HTTP = "1.10.19"
 JSON = "1.4.0"
 Minio = "0.2.2"
+OnlineStats = "1.7.3"
 Preferences = "1.5.1"
 ProgressMeter = "1.11.0"
 Proj = "1.9.0"
+SpectralIndices = "0.2.15"
 Statistics = "1.11.1"
 StringViews = "1.3.7"
-YAXArrays = "0.7.0"
+YAXArrays = "0.7.2"
 Zarr = "0.9.5"
 ZipArchives = "2.6.0"
 julia = "1.10"

src/ARCEMEAnalysis.jl

@@ -1,15 +1,29 @@
 module ARCEMEAnalysis
 using Zarr: S3Store, Zarr
 using Minio: MinioConfig
-using YAXArrays: open_dataset, ⊘, xmap, XOutput, YAXArray, YAXArrays, Dataset
+using YAXArrays: open_dataset, ⊘, xmap, XOutput, YAXArray, YAXArrays, Dataset, setchunks, compute_to_zarr, savedataset
 import DimensionalData as DD
 using Colors: RGB
 using Dates: DateTime, Year, Date
 import CSV
 using Statistics: mean
 using DataStructures: SortedDict, counter
 using ProgressMeter: @showprogress
+using SpectralIndices: compute_index, SpectralIndices
+import SpectralIndices as SI
+
+function __init__()
+    #Extend SpectralIndices band definitions with S1 bands
+    s1vv = Dict{String,SI.PlatformBand}("sentinel1" => SI.PlatformBand("sentinel1", "VV", "Vertical-Vertical", 5.55e7, 1e5))
+    s1vh = Dict{String,SI.PlatformBand}("sentinel1" => SI.PlatformBand("sentinel1", "VH", "Vertical-Horizontal", 5.55e7, 1e5))
+
+    SI.bands["VV"] = SI.Band("VV", "Vertical-Vertical", "vv", 5.54e7, 5.56e7, s1vv)
+    SI.bands["VH"] = SI.Band("VH", "Vertical-Horizontal", "vh", 5.54e7, 5.56e7, s1vh)
+end
+
+
 include("download.jl")
+include("s1_helpers.jl")
 
 const arceme_classes = SortedDict(
   0   => "No data",
@@ -25,11 +39,12 @@ const arceme_classes = SortedDict(
   95  => "Mangroves",
   100 => "Moss and lichen",
 )
+lckeymap(k) = ifelse(k > 90, (k + 10), k) ÷ 10 + 1
 
 export arceme_cubename, arceme_open, arceme_starttime, arceme_endtime, arceme_eventdate,
     arceme_coordinates, arceme_ndvi, arceme_rgb, arceme_eventlist, arceme_eventpairs, 
     arceme_classes, arceme_landcover, arceme_optical_band_fingerprints, arceme_radar_fingerprints,
-    time_aggregate_fingerprint, arceme_validpairs
+    time_aggregate_fingerprint, arceme_validpairs, arceme_spectral, arceme_kndvi, arceme_radar_db
 
 """
     _arceme_cubenames(;batch="6")
@@ -38,7 +53,7 @@ List all available data cube names in the specified batch stored in the ARCEME S
 """
 function _arceme_cubenames(; batch="ARCEME-DC-6")
     cubenames = if local_cubepath === nothing
-        store = S3Store("$(batch)/", MinioConfig("https://s3.waw3-2.cloudferro.com/swift/v1"))
+        store = S3Store("$(batch)/", MinioConfig(httpstore))
 
         resp = Zarr.cloud_list_objects(store, batch)
         map(split(String(resp), "\n")) do p
@@ -113,22 +128,23 @@ arceme_eventpairs() =
     Iterators.partition(sort(arceme_eventlist(),by=i->(i.dhp_label,i.source)),2) |> collect
 
 """
-    arceme_validpairs(;batch="ARCEME-DC-6", store="https://s3.waw3-2.cloudferro.com/swift/v1")
+    arceme_validpairs(;batch="ARCEME-DC-6")
 
-Get valid pairs of ARCEME events from the data store. Currently only working over HTTP(S)
+Get valid pairs of ARCEME events from the local path (if set with arceme_set_localpath) or the http data store.
+Default httpstore is "https://s3.waw3-2.cloudferro.com/swift/v1". It can be reset with arceme_set_httpstore.
 """
-function arceme_validpairs(;batch="ARCEME-DC-6", store="https://s3.waw3-2.cloudferro.com/swift/v1")
+function arceme_validpairs(;batch="ARCEME-DC-6")
     allpairs = arceme_eventpairs()
     validpairs = if local_cubepath === nothing
-        map(x -> all([Zarr.is_zgroup(Zarr.HTTPStore("$store/$batch/$(arceme_cubename(i))"), "") for i in x]), allpairs)
+        map(x -> all([Zarr.is_zgroup(Zarr.HTTPStore("$httpstore/$batch/$(arceme_cubename(i))"), "") for i in x]), allpairs)
     else
         map(x -> all([isfile(joinpath(local_cubepath, batch, string(arceme_cubename(i), ".zip"))) for i in x]), allpairs)
     end
     allpairs[validpairs]
 end
 
-function arceme_open(event::Event)
-    arceme_open(arceme_cubename(event))
+function arceme_open(event::Event; batch="ARCEME-DC-6")
+    arceme_open(arceme_cubename(event); batch=batch)
 end
 
 """
@@ -143,19 +159,29 @@ function arceme_landcover(ds)
         (key=k, class=v, count=count, fraction=count/1000000)
     end
 end
-arceme_landcover(ev::Event) = arceme_landcover(arceme_open(ev))
+arceme_landcover(ev::Event; batch="ARCEME-DC-6") = arceme_landcover(arceme_open(ev, batch=batch))
 
 
 """ 
     arceme_open(cubename; batch="ARCEME-DC-6")
 
-Open the specified ARCEME data cube from the S3 bucket.
+Open the specified ARCEME data cube from the local path (if with arceme_set_localpath) or 
+    the httpstore (default "https://s3.waw3-2.cloudferro.com/swift/v1", reset with arceme_set_httpstore).
 """
 function arceme_open(cubename; batch="ARCEME-DC-6")
     if local_cubepath === nothing
-        open_dataset("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/$cubename", force_datetime=true)
+        open_dataset("$httpstore/$batch/$cubename", force_datetime=true)
     else
-        open_dataset(joinpath(local_cubepath, batch, string(cubename, ".zip")))
+        main_ds = open_dataset(joinpath(local_cubepath, batch, string(cubename, ".zip")))
+        if isfile(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
+            index_ds = open_dataset(joinpath(local_cubepath, "$batch-INDICES", string(cubename, ".zip")))
+            for (k, v) in (index_ds.cubes)
+                main_ds.cubes[k] = v
+            end
+            return main_ds
+        else
+            return main_ds
+        end
     end
 end
 
@@ -205,15 +231,47 @@ Compute the NDVI (Normalized Difference Vegetation Index) for the ARCEME data cu
 """
 function arceme_ndvi(ds)
     ndvi = broadcast(ds.B04, ds.B08, ds.cloud_mask, ds.SCL) do b4, b8, cl, scl
-        (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11))) && return NaN
-        fb4 = b4 / typemax(Int16) # not necessary, right?
-        fb8 = b8 / typemax(Int16)
+        _is_cloud(cl, scl) && return NaN
+        fb4 = boa(b4)
+        fb8 = boa(b8)
         (fb8 - fb4) / (fb8 + fb4)
     end
     ds.cubes[:ndvi] = ndvi
     ds
 end
 
+"""
+    arceme_spectral(ds, indices::Vector{String})
+
+Compute the listed indices using SpectralIndices.jl. Not working.
+"""
+function arceme_spectral(ds, indices::Vector{String}; platform="sentinel2") 
+    tmp = if platform == "sentinel2" || platform == "sentinel2a" || platform == "sentinel2b"
+        pl = platform == "sentinel2" ? "sentinel2a" : platform
+        _compute_indices(ds, indices, pl)
+    elseif platform=="sentinel1"
+        _compute_indices(ds, indices, platform)
+    else
+        error("platform $platform is not supported")
+    end
+    foreach(pairs(tmp)) do (k,v)
+        ds.cubes[k] = v
+    end
+    ds
+end
+
+_is_cloud(cl, scl) = (cl > 0 || (scl in (1, 3, 7, 8, 9, 10, 11)))
+
+
+
+"""
+    boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000)
+
+Compute the radiometric offsets for Sentinel 2 bands to get values at bottom of atmosphere.
+The default values are valid for data processed from baseline 04.00 (January 2022) onwards.
+The entire Sentinel-2 archive in CDSE has been reprocessed and is now available in baseline 05.xx, with consistent offset.
+"""
+boa(band; BOA_ADD_OFFSET = -1000, QUANTIFICATION_VALUE = 10000) = (band + BOA_ADD_OFFSET) / QUANTIFICATION_VALUE
 
 """
     arceme_rgb(ds)
@@ -222,10 +280,10 @@ Compute the RGB composite for the ARCEME data cube dataset `ds`.
 """
 arceme_rgb(ds) =
     broadcast(ds.B02, ds.B03, ds.B04) do b, g, r
-        m = typemax(Int16)
+        # m = typemax(Int16)
         # RGB(r / m * 4, g / m * 4, b / m * 4)
         # RGB(min(1.0,r / m *4) , min(1.0,g / m *4) , min(1.0,b / m *4))
-        RGB(r / m , g / m , b / m )
+        RGB(clamp(boa(r), 0, 1), clamp(boa(g), 0, 1), clamp(boa(b), 0, 1))
 end
 
 """
@@ -336,7 +394,40 @@ function time_aggregate_fingerprint(allbands, eventdate, banddim, timeaxis)
     YAXArray((allbands.lc, DD.Ti((fingerprint_timesteps .- 0.5) ./ step_per_year .* 12), banddim), cat(res..., dims=2))
 end
 
+"""
+For every valid event pair creates a and stores data cubes of a list of precomputed vegetation indices. For kNDVI 
+a shared sigma parameter per land cover class is computed.
+"""
+function arceme_create_indexcubes(; indices_s1=["DpRVIVV"], indices_s2=["NDVI", "NDWI", "EVI2", "NIRv", "NDMI", "NSDSI3", "WDRVI"])
+
+    for ev in arceme_validpairs()
+
+        ds_pair = arceme_open.(ev)
+
+        foreach(ds_pair) do ds
+            arceme_spectral(ds, indices_s1, platform="sentinel1")
+            arceme_spectral(ds, indices_s2, platform="sentinel2")
+            arceme_radar_db(ds)
+        end
+        ARCEMEAnalysis.arceme_kndvi_pair(ds_pair...)
+
+        fields_to_save = [indices_s1; indices_s2]
+
+        foreach(ds_pair, ev) do ds, event
+            output_base = "$local_cubepath/ARCEME-DC-6-INDICES"
+            name = arceme_cubename(event)
+            indexcube = setchunks(ds[fields_to_save], (500, 500, 25))
+            compute_to_zarr(indexcube, joinpath(output_base, name), custom_loopranges=(500, 500, 25), overwrite=true)
+            cube2 = setchunks(ds[["vv_db", "vh_db", "kNDVI"]], (500, 500, 25))
+            savedataset(cube2, path=joinpath(output_base, name), append=true)
+            run(Cmd(`zip -0 -r ../$(name).zip .`, dir=joinpath(output_base, name)))
+            rm(joinpath(output_base, name), recursive=true)
+        end
+    end
+end
+
 
 
 include("spatialdebias.jl")
+include("spectral_helpers.jl")
 end #module

src/download.jl

@@ -3,7 +3,22 @@ using ProgressMeter
 using HTTP, JSON
 using ZipArchives
 using Preferences: @set_preferences!, @load_preference
-export arceme_set_localpath
+export arceme_set_localpath, arceme_set_httpstore
+
+"""
+arceme_set_httpstore(httpstore)
+
+Sets the https store from which to retrieve the cubes.
+If no local preference is set, the default location is used: https://s3.waw3-2.cloudferro.com/swift/v1
+
+Setting this preference needs recompilation.
+"""
+function arceme_set_httpstore(store)
+    @set_preferences!("arceme_httpstore" => store)
+    @warn "Preferences changed! Restart Julia for this change to take effect."
+end
+
+const httpstore = @load_preference("arceme_httpstore", "https://s3.waw3-2.cloudferro.com/swift/v1")
 
 """
     arceme_set_localpath(path)
@@ -12,6 +27,7 @@ Sets the path to a local directory where downloaded copies of the arceme cubes a
 """
 function arceme_set_localpath(path)
     @set_preferences!("arceme_localpath" => path)
+    @warn "Preferences changed! Restart Julia for this change to take effect."
 end
 
 const local_cubepath = @load_preference("arceme_localpath")
@@ -20,17 +36,17 @@ function arceme_download_batch(batch="ARCEME-DC-6")
     if local_cubepath === nothing
         error("You need to set the local cube path first. Please run `arceme_localpath(path)` first.")
     end
-    aresp = HTTP.get("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/",query=Dict("format"=>"json","delimiter"=>"/"))
+    aresp = HTTP.get("$httpstore/$batch/",query=Dict("format"=>"json","delimiter"=>"/"))
     allarrays = map(i->strip(i["subdir"],'/'),JSON.parse(aresp.body))
 
     @showprogress for current_cube in allarrays
-        lresp = HTTP.get("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/",query=Dict("prefix"=>current_cube))
+        lresp = HTTP.get("$httpstore/$batch/",query=Dict("prefix"=>current_cube))
         files_in_cube = split(StringView(lresp.body),"\n")
 
         outfilename = joinpath(local_cubepath,batch,string(current_cube,".zip"))
         ZipWriter(outfilename) do w
             for f in files_in_cube
-                resp = HTTP.get("https://s3.waw3-2.cloudferro.com/swift/v1/$batch/$f");
+                resp = HTTP.get("$httpstore/$batch/$f");
                 f2 = joinpath(splitpath(f)[2:end]...)
                 zip_newfile(w, f2)
                 write(w,resp.body)

src/s1_helpers.jl

@@ -0,0 +1,45 @@
+import Dates: Millisecond
+import YAXArrays: YAXArray
+"""
+    arceme_s1_position(ds)
+
+Analyses the sentinel 1 time stamps of the datasets and groups them into series of time stamps whose difference
+is always a multiple of 6 days, so we can assume that all images with the same group tag are retrieved from the 
+same position. Store the groups in the dataset with the name `s1_postion`.
+"""
+function arceme_s1_position(ds)
+    ts = ds.time_sentinel_1_rtc.val
+    sixdays = 6*24*60*60*1000
+    group_offsets = [0]
+    groups = [[1]]
+    for its in 2:length(ts)
+        groupfound = false
+        for igroup in eachindex(groups)
+            fac = (Millisecond(ts[its]-ts[1]).value-group_offsets[igroup])/sixdays
+            if abs(round(fac)-fac) < 1e-6
+                push!(groups[igroup],its)
+                groupfound=true
+                break
+            end
+        end
+        if !groupfound
+            newoffset = Dates.Millisecond(ts[its]-ts[1]).value-floor(Int,fac)*sixdays
+            push!(group_offsets,newoffset)
+            push!(groups,[its])
+        end
+    end
+    positions = zeros(Int,length(ts))
+    for i in 1:length(groups)
+        positions[groups[i]] .= i
+    end
+    positions
+    ds.cubes[:s1_position] = YAXArray((ds.time_sentinel_1_rtc,),positions)
+    nothing
+end
+
+_to_db(x) = 10.0 * log10(x)
+
+function arceme_radar_db(ds)
+    ds.cubes[:vv_db] = _to_db.(ds.vv)
+    ds.cubes[:vh_db] = _to_db.(ds.vh)
+end