Skip nan time steps and classify the rqa trend results into UInt8 (#101)

* Check for nan time steps and skip them

* Add sentence about scope of the project

* 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.

* Make lower and upperbound kwargs

* Fix tests

* Fix tests

Author: felixcremer

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"]

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).

src/rqatrend.jl

@@ -8,8 +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...)
-    @show outpath
-    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
@@ -31,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")
@@ -48,10 +47,24 @@ 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
+    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,)))
+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
@@ -117,7 +130,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

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