From c14ffdc0513377f8f0f804a4e19e9051a855736a Mon Sep 17 00:00:00 2001
From: Felix Cremer <felix.cremer@bgc-jena.mpg.de>
Date: Thu, 9 Jan 2025 13:54:49 +0100
Subject: [PATCH] Add main function

This is currently just a dump and might not work yet.
---
 Project.toml            |  2 +-
 src/RQADeforestation.jl | 61 ++++++++++++++++++++++++++++++++++++++++-
 src/auxil.jl            |  9 +++---
 3 files changed, 66 insertions(+), 6 deletions(-)

diff --git a/Project.toml b/Project.toml
index 00a2989..43221bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,7 +38,7 @@ Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"
 
 [compat]
 julia = "1.10"
-YAXArrayBase = "=0.6"
+YAXArrayBase = "0.6, 0.7"
 Glob = "1"
 NetCDF="0.12"
 Zarr = "0.9"
diff --git a/src/RQADeforestation.jl b/src/RQADeforestation.jl
index 9153110..cbb23e4 100644
--- a/src/RQADeforestation.jl
+++ b/src/RQADeforestation.jl
@@ -1,5 +1,5 @@
 module RQADeforestation
-__precompile__(false)
+#__precompile__(false)
 using Dates
 using ArchGDAL: ArchGDAL as AG
 using Glob
@@ -13,4 +13,63 @@ export gdalcube
 include("auxil.jl")
 include("analysis.jl")
 include("timestats.jl")
+
+
+function main(tiles = ["E048N018T3"]; pol="VH", orbit="*", thresh=3.0)
+    indir = "/eodc/products/eodc.eu/S1_CSAR_IWGRDH/SIG0/"
+    continent = "EU"
+    folders = ["V01R01","V0M2R4", "V1M0R1", "V1M1R1", "V1M1R2"]
+    corruptedfiles = "corrupted_tiles.txt"
+    # TODO save the corrupt files to a txt for investigation
+    for tilefolder in tiles
+
+        filenamelist = [glob("$(sub)/*$(continent)*20M/$(tilefolder)/*$(pol)_$(orbit)*.tif", indir) for sub in folders]
+
+        allfilenames = collect(Iterators.flatten(filenamelist))
+
+
+        relorbits = unique([split(basename(x), "_")[5][2:end] for x in allfilenames])
+        @show relorbits
+        for relorbit in relorbits
+            for y in [2018,2019,2020,2021,2022, 2023]
+
+            filenames = allfilenames[findall(contains("$(relorbit)_E"), allfilenames)]
+            @time cube = gdalcube(filenames)
+
+                path = joinpath(YAXDefaults.workdir[], "$(tilefolder)_rqatrend_$(pol)_$(relorbit)_thresh_$(thresh)_year_$(y)")
+                @show path
+                ispath(path*".done") && continue
+                ispath(path*"_zerotimesteps.done") && continue
+
+                tcube = cube[Time=Date(y-1, 7,1)..Date(y+1,7,1)]
+                @show size(cube)
+                @show size(tcube)
+                if size(tcube, Ti) == 0
+                    touch(path*"_zerotimesteps.done")
+                    continue
+                end
+            try
+                @time rqatrend(tcube; thresh, outpath=path * ".zarr", overwrite=true)
+            catch e 
+                
+                if e.captured.ex isa ArchGDAL.GDAL.GDALError
+                    println("Found GDALError:")
+                    println(e.captured.ex.msg)
+                   continue
+                else 
+                    rethrow(e)
+                end
+            end
+                #=@everywhere begin
+                    fname = "$(VERSION)_$(getpid())_$(time_ns()).heapsnapshot"
+                    Profile.take_heap_snapshot(fname;streaming=true)
+                end
+                =#
+
+                touch(path * ".done")
+            end
+        end
+    end
+end
+
 end
\ No newline at end of file
diff --git a/src/auxil.jl b/src/auxil.jl
index 0bed2c3..099f53b 100644
--- a/src/auxil.jl
+++ b/src/auxil.jl
@@ -1,4 +1,4 @@
-using YAXArrayBase: GDALBand, GDALDataset, get_var_handle
+using YAXArrayBase:  backendlist, get_var_handle
 using DiskArrayTools
 using DiskArrays: DiskArrays, GridChunks
 using DimensionalData: DimensionalData as DD, X, Y
@@ -107,14 +107,15 @@ function gdalcube(filenames::AbstractVector{<:AbstractString})
 
     #datasets = AG.readraster.(sfiles)
     onefile = first(sfiles)
-    yax1 = GDALDataset(onefile)
+    gd = backendlist[:gdal]
+    yax1 = gd(onefile)
     onecube = Cube(onefile)
     #@show onecube.axes
     gdb = get_var_handle(yax1, "Gray")
 
-    @assert gdb isa GDALBand
+    #@assert gdb isa GDALBand
     all_gdbs = map(sfiles) do f
-        gd = BufferGDALBand{eltype(gdb)}(f, gdb.band, gdb.size, gdb.attrs, gdb.cs, Dict{Int,AG.IRasterBand}())
+        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)