Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 72 additions & 72 deletions src/auxil.jl
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -17,27 +17,27 @@
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)

Check warning on line 26 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L20-L26

Added lines #L20 - L26 were not covered by tests
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)

Check warning on line 33 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L30-L33

Added lines #L30 - L33 were not covered by tests
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)
for ia in eachindex(arrays_now)
DiskArrays.readblock!(arrays_now[ia], view(buf, :, :, ia), i1, i2)

Check warning on line 37 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L36-L37

Added lines #L36 - L37 were not covered by tests
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)))

Check warning on line 40 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L39-L40

Added lines #L39 - L40 were not covered by tests
end
end

Expand Down Expand Up @@ -76,8 +76,8 @@

@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

"""
Expand All @@ -103,7 +103,7 @@
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
Expand All @@ -123,12 +123,12 @@
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
else
group += 1
group += 1
groups[i] = group
end
end
Expand Down Expand Up @@ -166,61 +166,61 @@
#@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")

Check warning on line 191 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L191

Added line #L191 was not covered by tests
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)

Check warning on line 202 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L202

Added line #L202 was not covered by tests

onefile = first(sfiles)
gd = backendlist[:gdal]
yax1 = gd(onefile)
onecube = Cube(onefile)

Check warning on line 207 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L204-L207

Added lines #L204 - L207 were not covered by tests
#@show onecube.axes
gdb = get_var_handle(yax1, "Gray")

Check warning on line 209 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L209

Added line #L209 was not covered by tests

#@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}())

Check warning on line 213 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L212-L213

Added lines #L212 - L213 were not covered by tests
end
stacked_gdbs = diskstack(all_gdbs)
attrs = copy(gdb.attrs)

Check warning on line 216 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L215-L216

Added lines #L215 - L216 were not covered by tests
#attrs["add_offset"] = Float16(attrs["add_offset"])
if haskey(attrs, "scale_factor")
attrs["scale_factor"] = Float16(attrs["scale_factor"])

Check warning on line 219 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L218-L219

Added lines #L218 - L219 were not covered by tests
end
all_cfs = CFDiskArray(stacked_gdbs, attrs)
return YAXArray((onecube.axes..., taxis), all_cfs, onecube.properties)

Check warning on line 222 in src/auxil.jl

View check run for this annotation

Codecov / codecov/patch

src/auxil.jl#L221-L222

Added lines #L221 - L222 were not covered by tests
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()
Expand Down
14 changes: 7 additions & 7 deletions src/rqatrend.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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


Expand Down Expand Up @@ -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

Expand All @@ -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, γ)
Expand Down
Loading