masks features

Project.toml

@@ -6,8 +6,10 @@ version = "0.1.0"
 [deps]
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c"
 Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"
@@ -17,7 +19,8 @@ CondaPkg = "0.2"
 DataStructures = "0.18"
 ProgressMeter = "1.10"
 PythonCall = "0.9"
-YAXArrays = "0.5"
 Statistics = "1.10"
+YAXArrays = "0.5"
 Zarr = "0.9"
+SparseArrays = "1.10"
 julia = "1.10"

src/UnpackSinTiles.jl

@@ -4,16 +4,20 @@ module UnpackSinTiles
     using Zarr
     using Statistics
     using DataStructures
-    using PythonCall
     using ProgressMeter
+    using Dates
+    using SparseArrays
 
     hdf(f) = @pyconst(pyimport("pyhdf.SD").SD)(f)
     export hdf
 
     include("loadTile.jl")
     include("metadata.jl")
+    include("pixelOperations.jl")
+    include("modisTiles.jl")
 
     export openTile, loadTileVariable
     export getTileKeys, getTilePath
+    export burnTimeSpan, updateAfterBurn!, aggTile
 end
 

src/loadTile.jl

@@ -8,7 +8,7 @@ function getTilePath(tile, in_date, root_path)
     list_files = readdir(path_month)
     # find name file for the given tile
     index = findfirst(x -> occursin(tile, x), list_files)
-    if length(index)>0 # it should be only 1
+    if !isnothing(index) && length(index) == 1 # it should be only 1
         tile_full_name = list_files[index]
         full_path = joinpath(path_month, tile_full_name)
         return full_path
@@ -23,7 +23,11 @@ end
 """
 function openTile(tile, in_date, root_path)
     full_path = getTilePath(tile, in_date, root_path)
-    return hdf(full_path)
+    if !isnothing(full_path)
+        return hdf(full_path)
+    else
+        return nothing
+    end
 end
 
 """
@@ -40,19 +44,24 @@ end
 """
 function loadTileVariable(tile, in_date, root_path, variable; close_file = true)
     hdf_tile = openTile(tile, in_date, root_path)
-    hdf_data = hdf_tile.select(variable).get()
-    if close_file
-        hdf_tile.end() # close hdf tile
+    if !isnothing(hdf_tile)
+        hdf_data = hdf_tile.select(variable).get()
+        if close_file
+            hdf_tile.end() # close hdf tile
+        end
+        hdf_data_jl = pyconvert(Array, hdf_data)
+        return hdf_data_jl
+    else
+        return nothing
     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)
+        bb_dates = loadTileVariable(tile, months_r[i], root_path, variable) # `nothing` propagates 'till here!
         if !isnothing(bb_dates)
-            burnYear[i][:,:] = bb_dates
+            burnYear[i][:,:] .= bb_dates
         end
     end
     return burnYear

src/metadata.jl

@@ -1,4 +1,4 @@
-export parse_metadata
+export parse_metadata, getTileMetadata
 
 """
     parse_metadata(in_hdf; metadata_field = "StructMetadata.0")
@@ -43,4 +43,12 @@ function parse_metadata(in_hdf; metadata_field = "StructMetadata.0")
         end
     end
     return result
+end
+
+function getTileMetadata(tile, in_date, root_path)
+    o = openTile(tile, in_date, root_path)
+    meta = parse_metadata(o)
+    m = meta["GridStructure"]["GRID_1"]
+    o.end()
+    return m
 end

src/modisTiles.jl

@@ -0,0 +1,53 @@
+export getTileInterval
+export modisTiles
+
+# MODIS TILES
+function modisTiles()
+    hs = "h" .* lpad.(string.(0:35), 2, '0')
+    vs = "v" .* lpad.(string.(0:17), 2, '0')
+    return [h*v for v in vs, h in hs]
+end
+
+
+"""
+    getTileInterval(hv_tile::String, n_h::Int=36, n_v::Int=18, h_bound::UnitRange=1:1440, v_bound::UnitRange=1:720)
+
+Given a label in the format "hXXvYY", this function calculates the corresponding intervals
+in the specified ranges `h_bound` and `v_bound`.
+
+# Arguments:
+- `hv_tile::String`: The input label, expected to be in the format "hXXvYY".
+- `n_h::Int`: Number of intervals to divide the `h_bound` range (default is 36).
+- `n_v::Int`: Number of intervals to divide the `v_bound` range (default is 18).
+- `h_bound::UnitRange`: The full range of `h` (default is `1:1440`).
+- `v_bound::UnitRange`: The full range of `v` (default is `1:720`).
+
+# Returns:
+- A tuple of two ranges: the interval in `h_bound` and the interval in `v_bound`.
+
+"""
+function getTileInterval(hv_tile, n_h=36, n_v=18, h_bound=1:1440, v_bound=1:720)
+    # Extract the h and v indices from the label (assuming format "hXXvYY")
+    h_index = parse(Int, hv_tile[2:3])
+    v_index = parse(Int, hv_tile[5:6])
+
+    h_step = length(h_bound) ÷ n_h
+    v_step = length(v_bound) ÷ n_v
+
+    # Calculate the start and end of the interval for h and v
+    h_start = h_bound[1] + h_index * h_step
+    h_end = h_start + h_step - 1
+
+    v_start = v_bound[1] + v_index * v_step
+    v_end = v_start + v_step - 1
+
+    return CartesianIndices((v_start:v_end, h_start:h_end)) #h_start:h_end, v_start:v_end
+end
+
+# Example usage:
+# result = getTileInterval("h35v17")
+# result = getTileInterval("h01v00")
+# result = getTileInterval("h00v01")
+
+# a_test = rand(720, 1440)
+# a_test[result]

src/pixelOperations.jl

@@ -77,7 +77,7 @@ function updateAfterBurn!(beforeBurn, afterBurn; burnWindow = 30)
     return nothing
 end
 
-function averageWindows(sMatrix, n::Int)
+function averageWindows(sMatrix, n::Int; mask=nothing)
     rows, cols = size(sMatrix)
     # Check if the matrix dimensions are divisible by n
     if rows % n != 0 || cols % n != 0
@@ -86,21 +86,40 @@ function averageWindows(sMatrix, n::Int)
     # 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)
+    result = fill(NaN32, new_rows, new_cols) # ? NaN32 here and then missing ?
+    n_smask = if isnothing(mask)
+        new_rows * new_cols
+    end
     # 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 : μ 
+            if !isnothing(mask)
+                smask = @view mask[(i-1)*n+1:i*n, (j-1)*n+1:j*n]
+                window .*= smask
+                n_smask = sum(smask)
+            end
+            μ = sum(window) / n_smask # this is the mean per type, given by mask.
+            result[i, j] = iszero(μ) ? NaN32 : μ
         end
     end
     return result
 end
 
-function agg(μBurn, afterBurn, tsteps; res = 60)
+# TODO: include masks here! land, pfts.
+# MODIS: LAND_COVER TYPE GLOBAL 500m
+function aggTile(μBurn, afterBurn, tsteps; res = 60)
     @showprogress for t_step in 1:tsteps
         μBurn[:,:,t_step] = averageWindows(afterBurn[t_step], res)
     end
     return μBurn
+end
+
+
+function _compute_mask(m, to_vals)
+    return m .== to_vals
+end
+
+function compute_mask(m, to_vals::AbstractVector)
+    return m .∈ Ref(to_vals)
 end