From ea61cd5dc867377e65cddcf04194c75ddb7f4e3a Mon Sep 17 00:00:00 2001 From: Daniel Loos Date: Mon, 16 Jun 2025 09:09:43 +0200 Subject: [PATCH 1/3] Fix format --- src/auxil.jl | 138 +++++++++++++++++++++++++-------------------------- 1 file changed, 69 insertions(+), 69 deletions(-) diff --git a/src/auxil.jl b/src/auxil.jl index 2da4998..1ea8ea4 100644 --- a/src/auxil.jl +++ b/src/auxil.jl @@ -1,4 +1,4 @@ -using YAXArrayBase: backendlist, get_var_handle +using YAXArrayBase: backendlist, get_var_handle using DiskArrayTools using DiskArrays: DiskArrays, GridChunks using DiskArrayEngine: DiskArrayEngine as DAE @@ -17,27 +17,27 @@ struct LazyAggDiskArray{T,F,A} <: AbstractDiskArray{T,3} 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) +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) +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) + 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))) + vbuf = view(buf, :, :, 1:length(arrays_now)) + map!(a.f, view(aout, :, :, j), eachslice(vbuf, dims=(1, 2))) end end @@ -76,8 +76,8 @@ end @testitem "getdate" begin using Dates - @test RQADeforestation.getdate("sometext20200919T202020_somemoretext1234") == DateTime(2020,9,19, 20,20,20) - @test_throws Exception RQADeforestation.getdate("sometext") + @test RQADeforestation.getdate("sometext20200919T202020_somemoretext1234") == DateTime(2020, 9, 19, 20, 20, 20) + @test_throws Exception RQADeforestation.getdate("sometext") end """ @@ -128,7 +128,7 @@ function stackindices(times, timediff=200000) if period.value < timediff groups[i] = group else - group += 1 + group += 1 groups[i] = group end end @@ -166,61 +166,61 @@ function gdalcube(filenames::AbstractVector{<:AbstractString}, stackgroups=:dae) #@show sdates # Put the dates which are 200 seconds apart into groups 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 + 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 + 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] - 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) + taxis = DD.Ti(dates_grouped) + gcube = Cube(sfiles[1]) + return YAXArray((DD.dims(gcube)[1:2]..., taxis), aggdata, gcube.properties,) 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) - onecube = Cube(onefile) - #@show onecube.axes - gdb = get_var_handle(yax1, "Gray") - - #@assert gdb isa GDALBand - all_gdbs = map(sfiles) do f - BufferGDALBand{eltype(gdb)}(f, gdb.band, gdb.size, gdb.attrs, gdb.cs, Dict{Int,AG.IRasterBand}()) - end - stacked_gdbs = diskstack(all_gdbs) - attrs = copy(gdb.attrs) - #attrs["add_offset"] = Float16(attrs["add_offset"]) - if haskey(attrs, "scale_factor") - attrs["scale_factor"] = Float16(attrs["scale_factor"]) + #datasets = AG.readraster.(sfiles) + taxis = DD.Ti(sdates) + + onefile = first(sfiles) + gd = backendlist[:gdal] + yax1 = gd(onefile) + onecube = Cube(onefile) + #@show onecube.axes + gdb = get_var_handle(yax1, "Gray") + + #@assert gdb isa GDALBand + all_gdbs = map(sfiles) do f + BufferGDALBand{eltype(gdb)}(f, gdb.band, gdb.size, gdb.attrs, gdb.cs, Dict{Int,AG.IRasterBand}()) + end + stacked_gdbs = diskstack(all_gdbs) + attrs = copy(gdb.attrs) + #attrs["add_offset"] = Float16(attrs["add_offset"]) + if haskey(attrs, "scale_factor") + attrs["scale_factor"] = Float16(attrs["scale_factor"]) + end + all_cfs = CFDiskArray(stacked_gdbs, attrs) + return YAXArray((onecube.axes..., taxis), all_cfs, onecube.properties) 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() From 4d0a45c51d78dccf445e16f54cebe31fc7000b39 Mon Sep 17 00:00:00 2001 From: Daniel Loos Date: Mon, 16 Jun 2025 09:16:18 +0200 Subject: [PATCH 2/3] Fix format --- src/rqatrend.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index cc388cf..553ab58 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -57,7 +57,7 @@ function rqatrend_impl(data; thresh=2, border=10, theiler=1, metric=CheckedEucli # simplified implementation of https://stats.stackexchange.com/a/370175 and https://github.com/joshday/OnlineStats.jl/blob/b89a99679b13e3047ff9c93a03c303c357931832/src/stats/linreg.jl # x is the diagonal offset, y the percentage of local recurrence # we compute the slope of a simple linear regression with bias from x to y - xs = 1+theiler : length(data)-border + xs = 1+theiler:length(data)-border x_mean = mean(xs) xx_mean = sqmean_step1_range(xs) # mean(x*x for x in xs) @@ -70,18 +70,18 @@ function rqatrend_impl(data; thresh=2, border=10, theiler=1, metric=CheckedEucli n += 1.0 y = tau_rr(data, x; thresh, metric) y_mean = smooth(y_mean, y, inv(n)) - xy_mean = smooth(xy_mean, x*y, inv(n)) + xy_mean = smooth(xy_mean, x * y, inv(n)) end - A = SA_F64[ + A = SA_F64[ xx_mean x_mean - x_mean 1.0 + x_mean 1.0 ] b = SA_F64[xy_mean, y_mean] # OnlineStats uses `Symmetric(A) \ b`, however this does not work for StaticArrays # `cholesky(A) \ b` is recommended instead at discourse https://discourse.julialang.org/t/staticarrays-solve-symmetric-linear-system-seems-typeinstable/124634 # some timings show that there is no significant speedup when adding cholesky or doing plain static linear regression # hence leaving it out for now - return 1000.0*(A \ b)[1] # slope + return 1000.0 * (A\b)[1] # slope end @@ -123,7 +123,7 @@ function tau_rr(y, d; thresh=2, metric=CheckedEuclidean()) nominator += evaluate(metric, y[i], y[i+d]) <= _thresh denominator += 1 end - return nominator/denominator + return nominator / denominator end end @@ -135,7 +135,7 @@ function sqmean_step1_range(xs) end # assumes n is Int for optimal performance -sumofsquares(n) = n*(n+1)*(2*n+1)/6 +sumofsquares(n) = n * (n + 1) * (2 * n + 1) / 6 """ smooth(a, b, γ) From 8a44fd66b7e1a682612018463d0af3603961fd00 Mon Sep 17 00:00:00 2001 From: Daniel Loos Date: Mon, 16 Jun 2025 09:26:44 +0200 Subject: [PATCH 3/3] Fix index style --- src/auxil.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/auxil.jl b/src/auxil.jl index 1ea8ea4..9a6ac1c 100644 --- a/src/auxil.jl +++ b/src/auxil.jl @@ -33,7 +33,7 @@ function DiskArrays.readblock!(a::LazyAggDiskArray, aout, i::UnitRange{Int}...) 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) + for ia in eachindex(arrays_now) DiskArrays.readblock!(arrays_now[ia], view(buf, :, :, ia), i1, i2) end vbuf = view(buf, :, :, 1:length(arrays_now)) @@ -103,7 +103,7 @@ function grouptimes(times, timediff=200000) group = [1] groups = [group] - for i in 2:length(times) + for i in eachindex(times)[2:end] t = times[i] period = t - times[group[end]] if period.value < timediff @@ -123,7 +123,7 @@ function stackindices(times, timediff=200000) group = 1 groups[1] = group - for i in 2:length(times) + for i in eachindex(times)[2:end] period = times[i] - times[i-1] if period.value < timediff groups[i] = group