operations

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c"

src/UnpackSinTiles.jl

@@ -3,12 +3,14 @@ module UnpackSinTiles
     using YAXArrays
     using Zarr
     using Statistics
+    using DataStructures
     using PythonCall
 
     hdf(f) = @pyconst(pyimport("pyhdf.SD").SD)(f)
     export hdf
 
     include("loadTile.jl")
+    include("metadata.jl")
 
     export openTile, loadTileVariable
     export getTileKeys, getTilePath

src/loadTile.jl

@@ -45,4 +45,35 @@ function loadTileVariable(tile, in_date, root_path, variable; close_file = true)
         hdf_tile.end() # close hdf tile
     end
     return hdf_data
+end
+
+function loadTileBurntYear(months_r, tile, root_path; variable = "Burn Date")
+    burnYear = [spzeros(Int16, 2400,2400) for _ in 1:12]
+    for i in 1:12
+        bb_dates = loadTileVariable(tile, months_r[i], root_path, variable)
+        if !isnothing(bb_dates)
+            burnYear[i][:,:] = bb_dates
+        end
+    end
+    return burnYear
+end
+
+function dailySparse(year, tile, root_path; variable="Burn Date")
+    doy_month = monthIntervals(year)
+    months_r = initMonthsArray(year)
+    loadedYear = loadTileBurntYear(months_r, tile, root_path; variable)
+    doy_bool= map(doy_month) do (doy, indx_month)
+        ba = loadedYear[indx_month]
+        getPixelDayRaw(ba, doy)
+    end
+    return doy_bool
+end
+
+function burnTimeSpan(s_year, e_year, tile, root_path; variable = "Burn Date")
+    all_years = dailySparse(s_year, tile, root_path; variable)
+    @showprogress for year in s_year+1:e_year
+        n_year = dailySparse(year, tile, root_path; variable)
+        push!(all_years, n_year...)
+    end
+    return all_years
 end

src/metadata.jl

@@ -0,0 +1,46 @@
+export parse_metadata
+
+"""
+    parse_metadata(in_hdf; metadata_field = "StructMetadata.0")
+
+The `metadata_field` follows from the MODIS MCD64A1.061 specification product.
+"""
+function parse_metadata(in_hdf; metadata_field = "StructMetadata.0")
+    attr = in_hdf.attributes();
+    metadata_string = string.(attr[metadata_field])
+
+    result = OrderedDict{String, Any}()
+    stack = [result]
+
+    for line in eachsplit(metadata_string, '\n')
+        line = rstrip(line)
+        level = count(x -> x == '\t', line)
+        line = strip(line)
+
+        if line == "END"
+            break  # Stop parsing when "END" is encountered
+        elseif startswith(line, "GROUP=") || startswith(line, "OBJECT=")
+            key = split(line, '=', limit=2)[2]
+            new_dict = OrderedDict{String, Any}()
+            stack[end][key] = new_dict
+            push!(stack, new_dict)
+        elseif startswith(line, "END_GROUP=") || startswith(line, "END_OBJECT=")
+            pop!(stack)
+        elseif contains(line, "=")
+            key, value = split(line, '=', limit=2)
+            value = strip(value, ['"'])  # Remove quotation marks
+            # Try to parse as number if possible
+            try
+                if contains(value, ",")  # Handle tuples
+                    value = Tuple(parse(Float64, v) for v in split(value[2:end-1], ','))
+                else
+                    value = parse(Float64, value)
+                end
+            catch
+                # If parsing fails, keep it as a string
+            end
+            stack[end][key] = value
+        end
+    end
+    return result
+end

src/pixelOperations.jl

@@ -0,0 +1,106 @@
+function getPixelDay(ba, doy)
+    day_b = ba .== doy # get indices for day d
+    day_mask = ba .< 1 # get indices for water, unmmaped and unburned pixels
+    bmask = ba .* day_mask # compute mask
+    return ba .* day_b .+ bmask # combine burnt pixels for day d, plus mask
+end
+
+# this one is just for fire or not fire, Bool, true -> 1, false -> 0
+function getPixelDayRaw(ba, doy)
+    day_b = ba .== doy # get indices for day d
+    return day_b # ba .* day_b # ouput only 0 or 1 for that day
+end
+
+function monthIntervals(year)
+    days = Date(year):Day(1):Date("$(year)-12-31")
+    doy = dayofyear.(days)
+    meses = Dates.month.(days)
+    return Pair.(doy, meses)
+end
+
+function initMonthsArray(year)
+    return ["$(year).$(m).01" for m in ["0".* string.(1:9)..., "10", "11", "12"]]
+end
+
+function getPixelTime(bSpaceTime, i,j)
+    return [bSpaceTime[k][i,j] for k in eachindex(bSpaceTime)]
+end
+
+function cluster_entries(vec::Vector{T}, gap_length::Int) where T<:Number
+    result = Vector{Vector{Int}}()
+    n = length(vec)
+    current_group = Int[]
+    zero_count = 0
+    for i in 1:n
+        if vec[i] != zero(T)
+            if zero_count >= gap_length && !isempty(current_group)
+                push!(result, current_group)
+                current_group = Int[]
+            end
+            push!(current_group, i)
+            zero_count = 0
+        else
+            zero_count += 1
+        end
+    end
+    if !isempty(current_group)
+        push!(result, current_group)
+    end
+    return result
+end
+
+function normalize_groups(groups::Vector{Vector{Int}})
+    normPairs = map(groups) do group
+        group_length = length(group)
+        normalized_value = 1 ./ group_length # 1 .// group_length
+        Pair.(group, normalized_value)
+    end
+    return vcat(normPairs...)
+end
+
+function updateCounts!(afterBurn, nPairs, i, j)
+    for (k, v) in nPairs
+        afterBurn[k][i,j] = v
+    end
+    return nothing
+end
+
+function updateAfterBurn!(beforeBurn, afterBurn; burnWindow = 30)
+    @showprogress for i in 1:size(beforeBurn[1], 1)
+        for j in 1:size(beforeBurn[1], 2)
+            onePix = getPixelTime(beforeBurn, i, j)
+            gBurn = cluster_entries(onePix, burnWindow)
+            normedPairs = normalize_groups(gBurn)
+            updateCounts!(afterBurn, normedPairs, i, j)
+        end
+    end
+    return nothing
+end
+
+function averageWindows(sMatrix, n::Int)
+    rows, cols = size(sMatrix)
+    # Check if the matrix dimensions are divisible by n
+    if rows % n != 0 || cols % n != 0
+        error("Matrix dimensions must be divisible by n")
+    end
+    # Calculate the size of the output matrix
+    new_rows, new_cols = rows ÷ n, cols ÷ n
+    # Initialize the output matrix
+    result = fill(NaN32, new_rows, new_cols)
+    # Iterate over the windows and calculate averages
+    for i in 1:new_rows
+        for j in 1:new_cols
+            window = @view sMatrix[(i-1)*n+1:i*n, (j-1)*n+1:j*n]
+            μ = mean(window)
+            result[i, j] = iszero(μ) ? NaN32 : μ 
+        end
+    end
+    return result
+end
+
+function agg(μBurn, afterBurn, tsteps; res = 60)
+    @showprogress for t_step in 1:tsteps
+        μBurn[:,:,t_step] = averageWindows(afterBurn[t_step], res)
+    end
+    return μBurn
+end

src/saveTile.jl

@@ -0,0 +1,14 @@
+
+# maybe start first with a complete empty zarr store and then modify in-place instead of creating YAXArrays on the fly
+
+function appendTile(tile, agg_μBurn, xdim, ydim, tdim, zarrPath)
+    axlist = (
+        Dim{:x}(xdim),
+        Dim{:y}(ydim),
+        Dim{:time}(tdim),
+    )
+    yax = YAXArray(axlist, agg_μBurn)
+    ds = YAXArrays.Dataset(; (Symbol(tile) => yax,)...)
+    # append to file
+    savedataset(ds, path=zarrPath, driver=:zarr, append=true)
+end

test/load.jl

@@ -8,4 +8,8 @@ tile_path = getTilePath(tile, in_date, root_path)
 
 o = openTile(tile, in_date, root_path)
 k_vars = getTileKeys(o)
+
+meta = parse_metadata(o)
+meta["GridStructure"]["GRID_1"]
+
 o.end()