Aggregate different scenes with DAE

Project.toml

@@ -44,8 +44,8 @@ ArchGDAL = "0.10"
 ArgParse = "1"
 ConstructionBase = "1"
 DimensionalData = "0.29"
-DiskArrayEngine = "0.2"
-DiskArrayTools = "0.1.12"
+DiskArrayEngine = "0.1, 0.2"
+DiskArrayTools = "0.1"
 DiskArrays = "0.4"
 Distances = "0.10"
 Distributions = "0.25"

src/analysis.jl

@@ -1,7 +1,6 @@
-using RecurrenceAnalysis
+using RecurrenceAnalysis: RecurrenceAnalysis as RA
 #using MultipleTesting
 using Distances
-const RA = RecurrenceAnalysis
 """
 countvalid(xout, xin)
 

src/auxil.jl

@@ -1,10 +1,47 @@
 using YAXArrayBase:  backendlist, get_var_handle
 using DiskArrayTools
 using DiskArrays: DiskArrays, GridChunks
+using DiskArrayEngine: DiskArrayEngine as DAE
 using DimensionalData: DimensionalData as DD, X, Y
 using GeoFormatTypes
 #using PyramidScheme: PyramidScheme as PS
 using Rasters: Raster
+import DiskArrays: readblock!, IrregularChunks, AbstractDiskArray
+using StatsBase: rle
+using Statistics: mean
+
+#DiskArrays.readblock!(a::AbstractArray,aout,i::AbstractUnitRange...) = copyto!(aout,view(a,i...))
+
+struct LazyAggDiskArray{T,F,A} <: AbstractDiskArray{T,3}
+    f::F
+    arrays::A
+    inds::IrregularChunks
+    s::Tuple{Int,Int,Int}
+end
+function LazyAggDiskArray(f,arrays,groups)
+    allequal(size,arrays) || error("All Arrays must have same size")
+    allequal(eltype,arrays) || error("All Arrays must have same element type")
+    inds = IrregularChunks(;chunksizes=last(rle(groups)))
+    s = (size(first(arrays))...,length(inds))
+    T = Base.promote_op(f,Vector{eltype(first(arrays))})
+    LazyAggDiskArray{T,typeof(f),typeof(arrays)}(f,arrays,inds,s)
+end
+Base.size(a::LazyAggDiskArray) = a.s
+DiskArrays.haschunks(a::LazyAggDiskArray) = DiskArrays.haschunks(first(a.arrays))
+function DiskArrays.readblock!(a::LazyAggDiskArray,aout,i::UnitRange{Int}...)
+    i1,i2,itime = i
+    max_n_array = maximum(it->length(a.inds[it]),itime)
+    buf = zeros(eltype(first(a.arrays)),length(i1),length(i2),max_n_array)
+    for (j, it) in enumerate(itime)
+        arrays_now = a.arrays[a.inds[it]]
+        for ia in 1:length(arrays_now)
+            DiskArrays.readblock!(arrays_now[ia],view(buf,:,:,ia),i1,i2)
+        end
+        vbuf = view(buf,:,:,1:length(arrays_now))
+        map!(a.f,view(aout,:,:,j),eachslice(vbuf,dims=(1,2)))
+    end
+end
+
 
 struct BufferGDALBand{T} <: AG.DiskArrays.AbstractDiskArray{T,2}
     filename::String
@@ -56,7 +93,7 @@ where the difference between neighbouring elements are less than `timediff` mill
 This returns the indices of the subgroups as a vector of vectors.
 """
 function grouptimes(times, timediff=200000)
-    @assert sort(times) == times
+    @assert issorted(times)
     group = [1]
     groups = [group]
 
@@ -73,6 +110,24 @@ function grouptimes(times, timediff=200000)
     return groups
 end
 
+function stackindices(times, timediff=200000)
+    @assert issorted(times)
+    groups = zero(eachindex(times))
+    group = 1
+    groups[1] = group
+
+    for i in 2:length(times)
+        period = times[i] - times[i-1]
+        if period.value < timediff
+            groups[i] = group
+        else
+            group += 1 
+            groups[i] = group
+        end
+    end
+    return groups
+end
+
 #=
 function DiskArrays.readblock!(b::GDALBand, aout, r::AbstractUnitRange...)
    if !isa(aout,Matrix)
@@ -93,19 +148,52 @@ function DiskArrays.readblock!(b::GDALBand, aout, r::AbstractUnitRange...)
 end
 =#
 
-function gdalcube(filenames::AbstractVector{<:AbstractString})
+function gdalcube(filenames::AbstractVector{<:AbstractString}, stackgroups=:dae)
     dates = getdate.(filenames)
+    @show length(dates)
     # Sort the dates and files by DateTime
     p = sortperm(dates)
     sdates = dates[p]
     sfiles = filenames[p]
-    taxis = DD.Ti(sdates)
 
     #@show sdates
     # Put the dates which are 200 seconds apart into groups
-    #groupinds = grouptimes(sdates, 200000)
+    if stackgroups in [:dae, :lazyagg]
+    groupinds = grouptimes(sdates, 200000)
+    onefile = first(sfiles)
+    gd = backendlist[:gdal]
+    yax1 = gd(onefile)
+    #gdb = yax1["Gray"]
+    #onecube = Cube(onefile)
+    #@show onecube.axes
+    gdb = get_var_handle(yax1, "Gray")
+    gdbband = gdb.band
+    gdbsize = gdb.size
+    gdbattrs = gdb.attrs
+    gdbcs = gdb.cs
+    group_gdbs = map(sfiles) do f
+            BufferGDALBand{eltype(gdb)}(f, gdbband, gdbsize, gdbattrs, gdbcs, Dict{Int,AG.IRasterBand}())
+    end
+
+    cubelist = CFDiskArray.(group_gdbs, (gdbattrs,))
+    stackinds = stackindices(sdates)
+    aggdata = if stackgroups == :dae
+        gcube = diskstack(cubelist)
+        aggdata = DAE.aggregate_diskarray(gcube, mean ∘ skipmissing, (3=> stackinds,); strategy=:direct)
+    else
+        println("Construct lazy diskarray")
+        LazyAggDiskArray(mean ∘ skipmissing, cubelist, stackinds)
+    end
+#    data = DiskArrays.ConcatDiskArray(reshape(groupcubes, (1,1,length(groupcubes))))
+    dates_grouped = [sdates[group[begin]] for group in groupinds]
 
+    taxis = DD.Ti(dates_grouped)
+    gcube = Cube(sfiles[1])
+    return YAXArray((DD.dims(gcube)[1:2]..., taxis), aggdata, gcube.properties,)    
+else
     #datasets = AG.readraster.(sfiles)
+    taxis = DD.Ti(sdates)
+
     onefile = first(sfiles)
     gd = backendlist[:gdal]
     yax1 = gd(onefile)
@@ -125,6 +213,7 @@ function gdalcube(filenames::AbstractVector{<:AbstractString})
     end
     all_cfs = CFDiskArray(stacked_gdbs, attrs)
     return YAXArray((onecube.axes..., taxis), all_cfs, onecube.properties)
+end
     #datasetgroups = [datasets[group] for group in groupinds]
     #We have to save the vrts because the usage of nested vrts is not working as a rasterdataset
     #temp = tempdir()

src/main.jl

@@ -1,6 +1,7 @@
 using ArgParse
 using YAXArrays: YAXDefaults
 
+
 const argparsesettings = ArgParseSettings()
 
 ArgParse.parse_item(::Type{Date}, x::AbstractString) = Date(x)
@@ -28,7 +29,7 @@ ArgParse.parse_item(::Type{Date}, x::AbstractString) = Date(x)
 
     "--orbit", "-o"
     help = "One of: Orbit number, 'A' for ascending, 'D' for descending, '*' for all orbits"
-    default = "*"
+    default = "A"
 
     "--out-dir", "-d"
     help = "Path to output zarr dataset"
@@ -72,10 +73,12 @@ function main(;
     start_date::Date,
     end_date::Date,
     polarisation="VH",
-    orbit="*",
+    orbit="D",
     threshold=3.0,
-    folders=["V01R01", "V0M2R4", "V1M0R1", "V1M1R1", "V1M1R2"]
+    folders=["V1M0R1", "V1M1R1", "V1M1R2"],
+    stack=:dae
 )
+    in(orbit, ["A", "D"]) || error("Orbit needs to be either A or D")
     if isdir(indir) && isempty(indir)
         error("Input directory $indir must not be empty")
     end
@@ -98,15 +101,14 @@ function main(;
 
         filenamelist = [glob("$(sub)/*$(continent)*20M/$(tilefolder)/*$(polarisation)_$(orbit)*.tif", indir) for sub in folders]
         allfilenames = collect(Iterators.flatten(filenamelist))
-        @show allfilenames
 
         relorbits = unique([split(basename(x), "_")[5][2:end] for x in allfilenames])
         @show relorbits
         for relorbit in relorbits
             filenames = allfilenames[findall(contains("$(relorbit)_E"), allfilenames)]
-            @time cube = gdalcube(filenames)
+            @time cube = gdalcube(filenames, stack)
 
-            path = joinpath(YAXDefaults.workdir[], "$(tilefolder)_rqatrend_$(polarisation)_$(relorbit)_thresh_$(threshold)")
+            path = joinpath(YAXDefaults.workdir[], "$(tilefolder)_rqatrend_$(polarisation)_$(orbit)$(relorbit)_thresh_$(threshold)_year_$(y)")
             @show path
             ispath(path * ".done") && continue
             ispath(path * "_zerotimesteps.done") && continue