Core area metrics

src/LandscapeMetrics.jl

@@ -47,12 +47,6 @@ export largestpatchindex
 include("area_and_edge/radiusofgyration.jl")
 export radiusofgyration
 
-include("area_and_edge/totaledge.jl")
-export totaledge
-
-include("area_and_edge/edgedensity.jl")
-export edgedensity
-
 # Shape
 include("shape/paratio.jl")
 export paratio, perimeterarearatio
@@ -61,5 +55,27 @@ export shapeindex
 include("shape/fractal.jl")
 export fractaldimensionindex
 
+# Core Area
+include("core_area/corearea.jl")
+export core_area
+
+include("core_area/coreareaindex.jl")
+export core_area_index
+
+include("core_area/totalcorearea.jl")
+export total_core_area
+
+include("core_area/Ncore.jl")
+export count_core_areas
+
+include("core_area/coreareapercentage.jl")
+export core_area_percentage
+
+include("core_area/numberofdisjunctcorearea.jl")
+export number_of_disjunct_core_areas
+
+include("core_area/disjunctcoreareadensity.jl")
+export disjunct_core_area_density
+
 end # module LandscapeMetrics
 

src/core_area/Ncore.jl

@@ -0,0 +1,70 @@
+"""
+    count_core_areas(A::Matrix{Int}, depth::Int=0)
+
+Count the number of core areas in a binary matrix A after eroding it by the specified depth.
+
+"""
+
+
+function count_core_areas(l::Landscape, patch, depth::Int=0)
+
+    # Create a mask for the specified patch
+    core_mask = patches(l) .== patch
+
+    # Erode the core mask by the specified depth
+    nrows, ncols = size(core_mask)
+    for _ in 1:depth
+        up    = vcat(falses(1, ncols), core_mask[1:end-1, :])
+        down  = vcat(core_mask[2:end, :], falses(1, ncols))
+        left  = hcat(falses(nrows, 1), core_mask[:, 1:end-1])
+        right = hcat(core_mask[:, 2:end], falses(nrows, 1))
+        core_mask = core_mask .& up .& down .& left .& right
+    end
+
+    # Now we need to count the number of connected components in core_mask
+    labels = zeros(Int, nrows, ncols)
+    label = 0
+
+    """
+        flood_fill(i, j, lab)
+
+        Flood fill algorithm to label connected components in the core_mask.
+    """
+    function flood_fill(i, j, lab)
+        stack = [(i, j)]
+        while !isempty(stack)
+            x, y = pop!(stack)
+            if 1 <= x <= nrows && 1 <= y <= ncols && core_mask[x, y] && labels[x, y] == 0
+                labels[x, y] = lab
+                for (dx, dy) in ((-1,0), (1,0), (0,-1), (0,1))
+                    push!(stack, (x + dx, y + dy))
+                end
+            end
+        end
+    end
+
+    for i in 1:nrows
+        for j in 1:ncols
+            if core_mask[i, j] && labels[i, j] == 0
+                label += 1
+                flood_fill(i, j, label)
+            end
+        end
+    end
+    return label
+end
+
+@testitem "Core area is positive for multi-cell patch with depth 1" begin
+    A = [
+        1 1 1 0 0 0 0;
+        1 1 1 0 2 2 2;
+        1 1 1 1 1 1 2;
+        0 0 0 1 1 1 2; 
+        0 0 0 1 1 1 0;
+        0 0 0 0 0 0 0;
+        0 0 0 0 0 0 0
+    ]
+    L = Landscape(A)
+
+    @test count_core_areas(L, 2, 1) == 2
+end

src/core_area/corearea.jl

@@ -0,0 +1,70 @@
+"""
+    core_area(l::Landscape, patch, depth)
+
+Core area of a patch in the landscape, defined as the area remaining
+after removing all edge cells up to the specified depth.
+"""
+
+function core_area(l::Landscape, patch, depth)
+
+    # Making a mask with all the cells in the patch
+    patch_mask = patches(l) .== patch
+
+    # Making a copy of the patch mask that will be modified to be only the core cells
+    core_mask = copy(patch_mask)
+
+
+    for d in 1:depth
+
+        # Identifying the cells on the edge of the patch with the specified depth of edge value
+        boundary = falses(size(core_mask))
+        for i in 1:size(core_mask, 1)
+            for j in 1:size(core_mask, 2)
+                if core_mask[i, j]
+                    for (di, dj) in ((-1,0), (1,0), (0,-1), (0,1))
+                        ni, nj = i+di, j+dj
+                        if 1 <= ni <= size(core_mask,1) && 1 <= nj <= size(core_mask,2)
+                            if !core_mask[ni, nj]
+                                boundary[i, j] = true
+                            end
+                        else
+                            boundary[i, j] = true 
+                        end
+                    end
+                end
+            end
+        end
+        core_mask[boundary] .= false
+    end
+    # Calculating the core area
+    cell_area = side(l)^2
+    core_area = sum(core_mask) * cell_area
+    return core_area
+end
+
+@testitem "Core area is zero for single cell patch with depth 1" begin
+    A = [0 1 0; 0 0 0; 0 0 0]
+    L = Landscape(A)
+    @test core_area(L, 1, 1) == 0
+end
+
+@testitem "Core area is positive for multi-cell patch with depth 1" begin
+    A = [
+        0 1 1 1;
+        0 1 1 1;
+        0 1 1 1
+    ]
+    L = Landscape(A)
+    @test core_area(L, 1, 1) == 1.0
+end
+
+@testitem "Core area is zero for multi-cell patch with depth 2" begin
+    A = [
+        0 1 1 1 1;
+        0 1 1 1 1;
+        0 1 1 1 1
+    ]
+    L = Landscape(A)
+    @test core_area(L, 1, 1) == 2 
+end
+

src/core_area/coreareaindex.jl

@@ -0,0 +1,34 @@
+"""
+    core_area_index(l::Landscape, patch, depth)
+
+Core area index (%) of a given patch in the landscape at a specified depth.
+
+"""
+
+function core_area_index(l::Landscape, patch, depth)
+
+    return core_area(l, patch, depth) / area(l, patch) * 100.0
+
+end
+
+@testitem "We can compute the core area index for a patch" begin
+    A = [
+        0 1 1 1 1;
+        0 1 1 1 1;
+        0 1 1 1 1
+    ]
+    L = Landscape(A)
+    patches!(L)
+    cai = core_area_index(L, 1, 1)
+    @test cai == (2/12)*100.0
+end
+
+@testitem "Core area is zero for multi-cell patch with depth 2" begin
+    A = [
+        0 1 1 1 1;
+        0 1 1 1 1;
+        0 1 1 1 1
+    ]
+    L = Landscape(A)
+    @test core_area(L, 1, 1) == 2 
+end

src/core_area/coreareapercentage.jl

@@ -0,0 +1,37 @@
+"""
+    core_area_percentage(l::Landscape, depth)
+
+Percentage of the landscape that is core area, given a specified edge depth.
+"""
+
+function core_area_percentage(l::Landscape, depth)
+
+    # Get all patch ids
+    patch_ids = unique(patches(l))
+
+    # Compute total core area
+    total_core = 0.0
+
+    # Sum core area for each patch
+    for pid in patch_ids
+        if pid != 0
+            total_core += core_area(l, pid, depth)
+        end
+    end
+    return (total_core / totalarea(l)) * 100.0
+end
+
+@testitem "We can compute the core area percentage for the landscape" begin
+    A = [
+        0 0 0 0 0 0;
+        0 1 1 1 1 0;
+        0 1 1 1 1 0;
+        0 1 1 1 1 0;
+        0 0 0 0 0 0;
+        0 0 0 0 0 0
+    ]
+    L = Landscape(A)
+    patches!(L)
+    cap = core_area_percentage(L, 1)
+    @test cap == (2/36)*100.0
+end

src/core_area/disjunctcoreareadensity.jl

@@ -0,0 +1,66 @@
+"""
+    disjunct_core_area_density(l::Landscape, class_id, depth)
+
+Disjunct core area density for a given class in the landscape.
+"""
+function disjunct_core_area_density(l::Landscape, class_id, depth)
+
+    return number_of_disjunct_core_areas(l, class_id, depth) / totalarea(l)
+end
+
+@testitem "We can compute the disjunct core area density for a class" begin
+    A = [
+        0 1 1 1 1 1 0;
+        0 1 1 1 1 1 0;
+        0 1 1 1 1 1 0;
+        0 0 0 0 0 0 0;
+        0 0 0 0 1 1 1;
+        0 0 0 0 1 1 1;
+        1 1 1 0 0 0 0;
+        1 1 1 0 0 0 0;
+        1 1 1 0 0 0 0
+    ]
+    L = Landscape(A)
+    patches!(L)
+    nodca = disjunct_core_area_density(L, 1, 1)
+    @test nodca == 2 / 63.0
+end
+
+
+"""
+    number_of_disjunct_core_areas(l::landscape, depth)
+
+    Number of disjunct core areas in each patch of the landscape at a specified depth.
+
+"""
+
+function disjunct_core_area_density(l::Landscape, depth)
+    total_disjunct_core_areas = 0
+    unique_patches = unique(patches(l))
+    for patch in unique_patches
+        if patch != 0  # Assuming 0 is the background/no-data value
+            total_disjunct_core_areas += count_core_areas(l, patch, depth)
+        end
+    end
+    return total_disjunct_core_areas / totalarea(l)
+
+end
+
+@testitem "We can compute the disjunct core area density for a class" begin
+    A = [
+        0 1 1 1 1 1 0;
+        0 1 1 1 1 1 0;
+        0 1 1 1 1 1 0;
+        0 0 0 0 0 0 0;
+        0 0 0 0 1 1 1;
+        0 0 0 0 1 1 1;
+        1 1 1 0 2 2 2;
+        1 1 1 0 2 2 2;
+        1 1 1 0 2 2 2
+    ]
+    L = Landscape(A)
+    patches!(L)
+    nodca = disjunct_core_area_density(L, 1)
+    @test nodca == 4 / 63.0
+end
+