From 34561cfaebc19316abe28500c7d028caeadc8176 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Thu, 5 Jun 2025 16:36:20 +0200 Subject: [PATCH 1/6] Check for nan time steps and skip them --- src/rqatrend.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index 553ab58..913ea19 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -117,7 +117,7 @@ function tau_rr(y, d; thresh=2, metric=CheckedEuclidean()) nominator = 0 denominator = 0 @inbounds for i in 1:length(y)-d - if y[i] === missing || y[i+d] === missing + if y[i] === missing || y[i+d] === missing|| isnan(y[i]) || isnan(y[i+d]) continue end nominator += evaluate(metric, y[i], y[i+d]) <= _thresh From 48a05826938a4637ad31c31cfebaa37569b93dae Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 18 Aug 2025 14:51:11 +0200 Subject: [PATCH 2/6] Add sentence about scope of the project --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index c8f1ad8..2ccd02e 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,8 @@ This is the Julia package developed within the [FAIRSenDD project](https://github.com/EarthyScience/FAIRSenDD). +The aim of this project is to provide a forest change detection algorithm in Julia. + ## Citation F. Cremer, M. Urbazaev, J. Cortés, J. Truckenbrodt, C. Schmullius and C. Thiel, "Potential of Recurrence Metrics from Sentinel-1 Time Series for Deforestation Mapping," in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 13, pp. 5233-5240, 2020, doi: [10.1109/JSTARS.2020.3019333](https://dx.doi.org/10.1109/JSTARS.2020.3019333). From 251431bc7a001989c0e8b913c5174bb10638a817 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 8 Sep 2025 10:14:31 +0200 Subject: [PATCH 3/6] Add UInt8 classification This adds a classification into UInt8 buckets because we don't need the full accuracy of Float. This is a good compromise between flexibility and size of the output data. --- src/rqatrend.jl | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index 913ea19..65620d7 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -48,10 +48,23 @@ 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(pix_trend, pix, thresh=2) - pix_trend .= rqatrend_impl(pix; thresh) +function rqatrend(pix_trend, pix, thresh=2, lowerbound=-5., upperbound=-0.5) + pix_trend .= classify_rqatrend(rqatrend_impl(pix; thresh); lowerbound, upperbound) end +function classify_rqatrend(trend, lowerbound=Float32(-5.0), upperbound=Float32(-0.5)) + ctrend = clamp(trend, lowerbound, upperbound) + rlength = upperbound - lowerbound + return round(UInt8, 255-((ctrend - lowerbound) / rlength) * 255) +end + +@testitem "classify_rqatrend" begin + @test RQADeforestation.classify_rqatrend(-4.999) === UInt8(255) + @test RQADeforestation.classify_rqatrend(1) === UInt8(0) + @test RQADeforestation.classify_rqatrend(-0.52) === UInt8(1) + @test RQADeforestation.classify_rqatrend(-6) === UInt8(255) + @test isempty( AllocCheck.check_allocs(RQADeforestation.classify_rqatrend, (Float32, Float32, Float32))) +end function rqatrend_impl(data; thresh=2, border=10, theiler=1, metric=CheckedEuclidean()) # simplified implementation of https://stats.stackexchange.com/a/370175 and https://github.com/joshday/OnlineStats.jl/blob/b89a99679b13e3047ff9c93a03c303c357931832/src/stats/linreg.jl From 07ff66b23ce9f75f94da24b1f48d5cc606ebbbfb Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 8 Sep 2025 10:39:50 +0200 Subject: [PATCH 4/6] Make lower and upperbound kwargs --- src/rqatrend.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index 65620d7..54f5a41 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -8,7 +8,6 @@ using Distances Compute the RQA trend metric for the datacube `cube` with the epsilon threshold `thresh`. """ function rqatrend(cube; thresh=2, outpath=tempname() * ".zarr", overwrite=false, kwargs...) - @show outpath mapCube(rqatrend, cube, thresh; indims=InDims("Time"), outdims=OutDims(; outtype=Float32, path=outpath, overwrite, kwargs...)) end @@ -52,7 +51,7 @@ function rqatrend(pix_trend, pix, thresh=2, lowerbound=-5., upperbound=-0.5) pix_trend .= classify_rqatrend(rqatrend_impl(pix; thresh); lowerbound, upperbound) end -function classify_rqatrend(trend, lowerbound=Float32(-5.0), upperbound=Float32(-0.5)) +function classify_rqatrend(trend; lowerbound=Float32(-5.0), upperbound=Float32(-0.5)) ctrend = clamp(trend, lowerbound, upperbound) rlength = upperbound - lowerbound return round(UInt8, 255-((ctrend - lowerbound) / rlength) * 255) From d01d3d21bcef09a2c399d224ec544dfc68dec533 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 8 Sep 2025 11:29:49 +0200 Subject: [PATCH 5/6] Fix tests --- src/rqatrend.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/rqatrend.jl b/src/rqatrend.jl index 54f5a41..16a8dbb 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -8,7 +8,7 @@ using Distances Compute the RQA trend metric for the datacube `cube` with the epsilon threshold `thresh`. """ function rqatrend(cube; thresh=2, outpath=tempname() * ".zarr", overwrite=false, kwargs...) - mapCube(rqatrend, cube, thresh; indims=InDims("Time"), outdims=OutDims(; outtype=Float32, path=outpath, overwrite, kwargs...)) + mapCube(rqatrend, cube, thresh; indims=InDims("Time"), outdims=OutDims(; outtype=UInt8, path=outpath, overwrite, kwargs...)) end @testitem "rqatrend cube" begin @@ -30,8 +30,8 @@ end mock_trend = rqatrend(mock_cube; thresh=0.5) @test mock_trend.axes == (mock_cube.X, mock_cube.Y) - diff = abs(mean(mock_trend)) - @test diff < 0.5 + @test eltype(mock_trend) == Union{Missing, UInt8} + end """rqatrend(path::AbstractString; thresh=2, outpath=tempname()*".zarr") From db5a2abb9cb45fd33984c7358e3d60f66b3cb6ed Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 8 Sep 2025 14:36:56 +0200 Subject: [PATCH 6/6] Fix tests --- Project.toml | 6 ++++-- src/rqatrend.jl | 3 ++- test/testdata.jl | 3 +-- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index b9471f1..8f2abac 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a" GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" KML = "1284bf3a-1e3d-4f4e-a7a9-b9d235a28f35" +LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" @@ -61,6 +62,7 @@ GDAL = "1" GeoFormatTypes = "0.4" Glob = "1" KML = "0.2" +LazyArtifacts = "1.11.0" Libdl = "1.10" LinearAlgebra = "1.10" LoggingExtras = "1" @@ -82,7 +84,7 @@ TestItems = "1.0" TimeseriesSurrogates = "2" UnicodePlots = "3" YAXArrayBase = "0.6, 0.7" -YAXArrays = "0.5, 0.6" +YAXArrays = "0.5, 0.6, 0.7" Zarr = "0.9" julia = "1.11" @@ -91,13 +93,13 @@ AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StaticTools = "86c06d3c-3f03-46de-9781-57580aa96d0a" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" -Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" [targets] test = ["Test", "TestItemRunner", "Pkg", "Random", "AllocCheck", "BenchmarkTools", "Aqua", "Documenter", "StaticTools", "PythonCall", "Libdl"] diff --git a/src/rqatrend.jl b/src/rqatrend.jl index 16a8dbb..17967aa 100644 --- a/src/rqatrend.jl +++ b/src/rqatrend.jl @@ -58,11 +58,12 @@ function classify_rqatrend(trend; lowerbound=Float32(-5.0), upperbound=Float32(- end @testitem "classify_rqatrend" begin + import AllocCheck @test RQADeforestation.classify_rqatrend(-4.999) === UInt8(255) @test RQADeforestation.classify_rqatrend(1) === UInt8(0) @test RQADeforestation.classify_rqatrend(-0.52) === UInt8(1) @test RQADeforestation.classify_rqatrend(-6) === UInt8(255) - @test isempty( AllocCheck.check_allocs(RQADeforestation.classify_rqatrend, (Float32, Float32, Float32))) + @test isempty( AllocCheck.check_allocs(RQADeforestation.classify_rqatrend, (Float32,))) end function rqatrend_impl(data; thresh=2, border=10, theiler=1, metric=CheckedEuclidean()) diff --git a/test/testdata.jl b/test/testdata.jl index 2b548a1..cd2a9b5 100644 --- a/test/testdata.jl +++ b/test/testdata.jl @@ -1,6 +1,7 @@ @testitem "testdata main" begin import Pkg: Artifacts.@artifact_str + using LazyArtifacts testdatapath = artifact"rqatestdata/RQADeforestationTestData-2.0" testdir = tempname() @@ -24,8 +25,6 @@ a = open_dataset(outdir * "/E051N018T3_rqatrend_VH_D022_thresh_3.0.zarr").layer @test size(a) == (50, 74) - @test minimum(a) < 0 - @test maximum(a) > 0 end @testitem "testdata julia_main" begin