diff --git a/src/RQADeforestation.jl b/src/RQADeforestation.jl index c409457..9587b48 100644 --- a/src/RQADeforestation.jl +++ b/src/RQADeforestation.jl @@ -15,7 +15,6 @@ include("metrics.jl") include("auxil.jl") include("rqatrend.jl") include("analysis.jl") # TODO what is still needed from analysis now that rqatrend is in its own file? -include("timestats.jl") include("main.jl") end \ No newline at end of file diff --git a/src/analysis.jl b/src/analysis.jl index b1c8131..da2853d 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -35,83 +35,16 @@ countvalid(cube; path=tempname() * ".zarr") = mapCube(countvalid, cube; indims=I Y(range(1, 5, length=15)), ) mock_data = allowmissing(rand(30, 10, 15)) - mock_data[1:10,1,1] .= missing - mock_data[:, 2,1] .= missing - mock_data[[1,5,9], 2,2] .= missing + mock_data[1:10, 1, 1] .= missing + mock_data[:, 2, 1] .= missing + mock_data[[1, 5, 9], 2, 2] .= missing mock_props = Dict() mock_cube = YAXArray(mock_axes, mock_data, mock_props) mock_count = RQADeforestation.countvalid(mock_cube) @test mock_count.axes == (mock_cube.X, mock_cube.Y) - @test mock_count[1,1] == 20 - @test mock_count[1,2] == 30 - @test mock_count[2,2] == 27 - @test mock_count[2,1] == 0 -end - -""" -rqatrend(xout, xin, thresh) - -Compute the RQA trend metric for the non-missing time steps of xin, and save it to xout. -`thresh` specifies the epsilon threshold of the Recurrence Plot computation -""" -function rqatrend_recurrenceanalysis(pix_trend, pix, thresh=2) - ts = collect(skipmissing(pix)) - tau_pix = tau_recurrence(ts, thresh) - pix_trend .= RA._trend(tau_pix) -end - -function rqatrend_matrix(pix_trend, pix, thresh=2) - #replace!(pix, -9999 => missing) - ts = collect(skipmissing(pix)) - rm = RecurrenceMatrix(ts, thresh) - pix_trend .= RA.trend(rm) -end - -#= -""" - rqatrend_shuffle(cube; thresh=2, path=tempname() * ".zarr", numshuffle=300) -Compute the RQA trend metric for shuffled time series of the data cube `cube` with the epsilon threshold `thresh` for `numshuffle` tries and save it into `path`. -""" -function rqatrend_shuffle(cube; thresh=2, path=tempname() * ".zarr", numshuffle=300) - # This should be made a random shuffle - # TODO this looks completely broken - sg = surrogenerator(collect(eachindex(water[overlap])), BlockShuffle(7, shift=true)) -end -=# - - -""" - anti_diagonal_density(ts, thresh, metric) -Compute the average density of the diagonals perpendicular to the main diagonal for data series `ts`. -Uses the threshold `thresh` and `metric` for the computation of the similarities. -""" -function anti_diagonal_density(ts::AbstractVector, thresh, metric=Euclidean()) - n = length(ts) - ad_densities = zeros(2 * n - 3) - for col in 1:n - for row in 1:(col-1) - d = evaluate(metric, ts[col], ts[row]) - #@show row, col, d - ad_densities[col+row-2] += d <= thresh - end - end - half = div(n, 2) - maxdensities = collect(Iterators.flatten([(n, n) for n in 1:half-1])) - diagonallengths = [maxdensities..., half, reverse(maxdensities)...] - ad_densities ./ diagonallengths -end - -""" -Compute the forest masking thresholding and clustering of the rqadata in one step -""" -function inner_postprocessing(rqadata, forestmask; threshold=-1.28, clustersize=30) - @time rqamasked = rqadata .* forestmask - @time rqathresh = map(rqamasked) do x - if !ismissing(x) - x > threshold ? zero(Float32) : one(Float32) - else - x - end - end + @test mock_count[1, 1] == 20 + @test mock_count[1, 2] == 30 + @test mock_count[2, 2] == 27 + @test mock_count[2, 1] == 0 end \ No newline at end of file diff --git a/src/timestats.jl b/src/timestats.jl deleted file mode 100644 index 5591845..0000000 --- a/src/timestats.jl +++ /dev/null @@ -1,37 +0,0 @@ -using Statistics -using StatsBase - -function timestats(cube; path=tempname()) - - indims = InDims("Time") - funcs = ["Mean", "5th Quantile", "25th Quantile", "Median", "75th Quantile", "95th Quantile", - "Standard Deviation", - "Minimum", "Maximum", - "Skewness", "Kurtosis", "Median Absolute Deviation"] - - stataxis = CategoricalAxis("Stats", funcs) - od = OutDims(stataxis, path=path) - stats = mapCube(ctimestats!, cube, indims=indims, outdims=od) -end - -function ctimestats!(xout, xin) - x = collect(skipmissing(xin)) - ts = x[.!isnan.(x)] - # m = mean(ts) - # T = eltype(m) - #stats = Vector{T}(undef,12) - if isempty(ts) - xout .= NaN - else - xout[1] = mean(ts) - xout[2:6] .= quantile(ts, [0.05, 0.25, 0.5, 0.75, 0.95]) - xout[7] = std(ts) - xout[2] = minimum(ts) - xout[3] = maximum(ts) - # stats[10] = skewness(ts) - # stats[11] = kurtosis(ts) - # stats[12] = mad(ts, normalize=true) - end - #xout .=stats - nothing -end \ No newline at end of file diff --git a/test/Artifacts.toml b/test/Artifacts.toml index 9ba9b91..db39948 100644 --- a/test/Artifacts.toml +++ b/test/Artifacts.toml @@ -2,5 +2,5 @@ git-tree-sha1 = "1dc9d7aa975de4f5ecc652f34d6bbee706fcb7bd" [[rqatestdata.download]] - url = "https://github.com/meggart/RQADeforestationTestData/archive/refs/tags/v2.0.tar.gz" + url = "https://github.com/EarthyScience/RQADeforestationTestData/archive/refs/tags/v2.0.tar.gz" sha256 = "38a3e66a0505ed9d22c6e065627039089a421b0b143ac4c8d78aac45dd3dd497"