From b86a6ccda52929308588d247100eb317f51d4b30 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Wed, 19 Nov 2025 14:01:51 +0100 Subject: [PATCH 1/3] Add percentile range --- src/RQADeforestation.jl | 2 +- src/main.jl | 12 ++++++++++-- src/rqatrend.jl | 39 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 50 insertions(+), 3 deletions(-) diff --git a/src/RQADeforestation.jl b/src/RQADeforestation.jl index f2d4664..b37ffd8 100644 --- a/src/RQADeforestation.jl +++ b/src/RQADeforestation.jl @@ -10,7 +10,7 @@ using NetCDF using TestItems using DimensionalData: DimensionalData as DD -export gdalcube, rqatrend +export gdalcube, rqatrend, prange include("metrics.jl") include("auxil.jl") diff --git a/src/main.jl b/src/main.jl index 7a9cb03..b84f360 100644 --- a/src/main.jl +++ b/src/main.jl @@ -78,7 +78,8 @@ function main(; threshold=3.0, folders=["V1M0R1", "V1M1R1", "V1M1R2"], stack=:lazyagg, - delete_intermediate=false + delete_intermediate=false, + compute_prange=true, ) outdir = Path(outdir) in(orbit, ["A", "D"]) || error("Orbit needs to be either A or D") @@ -140,7 +141,14 @@ function main(; c = Cube(tmppath) @time "save to S3" savecube(setchunks(c, (15000,15000)), string(orbitoutpath)) rm(tmppath, recursive=true) - @show delete_intermediate + if compute_prange + prangepath = replace(orbitoutpath, "rqatrend"=>"prange", "thresh_3.0_"=>"") + tmppath = tempname() * ".zarr" + @time "prange" prange(tcube, outpath=tmppath, overwrite=true) + cprange = Cube(tmppath) + @time "save to S3" savecube(setchunks(c, (15000,15000)), string(prangepath)) + rm(tmppath, recursive=true) + end if delete_intermediate == false #PyramidScheme.buildpyramids(orbitoutpath) Zarr.consolidate_metadata(string(orbitoutpath)) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index b2dcb3b..41ae50a 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -170,3 +170,42 @@ Weighted average of `a` and `b` with weight `γ`. ``(1 - γ) * a + γ * b`` """ smooth(a, b, γ) = a + γ * (b - a) + +""" +prange(xout, xin) + +Compute the percentile range for the non-missing time steps of xin, and save it to xout. +`lowerpercentile` and `upperpercentile` specify the boundary of the percentile range. +These have to be between 0 and 1. +""" +function prange(xout, xin, lowpercentile=0.02, upperpercentile=0.98) + lowerlim, upperlim = quantile(xin, [lowpercentile, upperpercentile]) + xout .= upperlim - lowerlim +end + +function prange(cube; lowerpercentile=0.02, upperpercentile=0.98, outpath=tempname() * ".zarr", overwrite=false, kwargs...) + mapCube(prange, cube, lowerpercentile, upperpercentile; indims=InDims("Time"), outdims=OutDims(; outtype=Float32, path=outpath, fill_value=NaN, overwrite, kwargs...)) +end + +@testitem "prange cube" begin + using YAXArrays + using Dates + using DimensionalData: Ti, X, Y + using Statistics + import Random + Random.seed!(1234) + + mock_axes = ( + Ti(Date("2022-01-01"):Day(1):Date("2022-01-01")+ Day(100)), + X(range(1, 10, length=10)), + Y(range(1, 5, length=15)), + ) + s = size(mock_axes) + mock_data = reshape(1:prod(s), s) + mock_props = Dict() + mock_cube = YAXArray(mock_axes, mock_data, mock_props) + + mock_trend = prange(mock_cube) + @test mock_trend.axes == (mock_cube.X, mock_cube.Y) + @test mock_trend[1,1] == 96 +end \ No newline at end of file From aaf200979a1604e333afd4dc80d26f733037f100 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Thu, 20 Nov 2025 13:46:37 +0100 Subject: [PATCH 2/3] Fix replace on a path --- src/main.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.jl b/src/main.jl index b84f360..f3ddf2f 100644 --- a/src/main.jl +++ b/src/main.jl @@ -142,7 +142,7 @@ function main(; @time "save to S3" savecube(setchunks(c, (15000,15000)), string(orbitoutpath)) rm(tmppath, recursive=true) if compute_prange - prangepath = replace(orbitoutpath, "rqatrend"=>"prange", "thresh_3.0_"=>"") + prangepath = Path(replace(string(orbitoutpath), "rqatrend"=>"prange", "thresh_3.0_"=>"")) tmppath = tempname() * ".zarr" @time "prange" prange(tcube, outpath=tmppath, overwrite=true) cprange = Cube(tmppath) From 70986d01efae2f154fc56d64a1dfef4e563cbd27 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Fri, 28 Nov 2025 09:50:28 +0100 Subject: [PATCH 3/3] Filter nans and missings inside prange inner function --- src/rqatrend.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index 41ae50a..44cfb78 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -179,7 +179,9 @@ Compute the percentile range for the non-missing time steps of xin, and save it These have to be between 0 and 1. """ function prange(xout, xin, lowpercentile=0.02, upperpercentile=0.98) - lowerlim, upperlim = quantile(xin, [lowpercentile, upperpercentile]) + xinfiltered = filter(!ismissing, xin) + filter!(!isnan, xinfiltered) + lowerlim, upperlim = quantile(xinfiltered, [lowpercentile, upperpercentile]) xout .= upperlim - lowerlim end