From a9241f2402674528fb933c13583f4a006cbbb04b Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Sun, 1 Dec 2024 13:00:39 +0100
Subject: [PATCH 01/11] start Rasters.jl compatibility

---
 Project.toml   |  2 ++
 src/grid.jl    | 23 ++++++++++++++++++-----
 src/gridrsp.jl | 26 +++++++++++++++++---------
 src/utils.jl   | 12 +++++-------
 4 files changed, 42 insertions(+), 21 deletions(-)

diff --git a/Project.toml b/Project.toml
index 872d0e3..274de50 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
+Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -21,5 +22,6 @@ Graphs = "1"
 LaTeXStrings = "1.1"
 Plots = "1.4"
 ProgressLogging = "0.1"
+Rasters = "0.12"
 SimpleWeightedGraphs = "1.1"
 julia = "1.9"
diff --git a/src/grid.jl b/src/grid.jl
index e79e676..f31863e 100644
--- a/src/grid.jl
+++ b/src/grid.jl
@@ -26,6 +26,7 @@ struct Grid
     id_to_grid_coordinate_list::Vector{CartesianIndex{2}}
     source_qualities::Matrix{Float64}
     target_qualities::AbstractMatrix{Float64}
+    dims::Tuple
 end
 
 """
@@ -47,8 +48,8 @@ affinity and cost matrices will be pruned to exclude unreachable nodes.
 function Grid(nrows::Integer,
               ncols::Integer;
               affinities=nothing,
-              qualities::Matrix=ones(nrows, ncols),
-              source_qualities::Matrix=qualities,
+              qualities::AbstractMatrix=ones(nrows, ncols),
+              source_qualities::AbstractMatrix=qualities,
               target_qualities::AbstractMatrix=qualities,
               costs::Union{Transformation,SparseMatrixCSC{Float64,Int}}=MinusLog(),
               prune=true)
@@ -103,6 +104,7 @@ function Grid(nrows::Integer,
         id_to_grid_coordinate_list,
         _source_qualities,
         _target_qualities,
+        dims(source_qualities),
     )
 
     if prune
@@ -113,6 +115,7 @@ function Grid(nrows::Integer,
 end
 
 Base.size(g::Grid) = (g.nrows, g.ncols)
+DimensionalData.dims(g::Grid) = g.dims
 
 function Base.show(io::IO, ::MIME"text/plain", g::Grid)
     print(io, summary(g), " of size ", g.nrows, "x", g.ncols)
@@ -155,7 +158,7 @@ function plot_values(g::Grid, values::Vector; kwargs...)
     for (i,v) in enumerate(values)
         canvas[g.id_to_grid_coordinate_list[i]] = v
     end
-    heatmap(canvas, yflip=true, axis=nothing, border=:none, aspect_ratio=:equal; kwargs...)
+    _heatmap(canvas, g; kwargs...)
 end
 
 function plot_outdegrees(g::Grid; kwargs...)
@@ -164,7 +167,7 @@ function plot_outdegrees(g::Grid; kwargs...)
     for (i,v) in enumerate(values)
         canvas[g.id_to_grid_coordinate_list[i]] = v
     end
-    heatmap(canvas, yflip=true, axis=nothing, border=:none; kwargs...)
+    _heatmap(canvas, g; kwargs...)
 end
 
 function plot_indegrees(g::Grid; kwargs...)
@@ -173,7 +176,17 @@ function plot_indegrees(g::Grid; kwargs...)
     for (i,v) in enumerate(values)
         canvas[g.id_to_grid_coordinate_list[i]] = v
     end
-    heatmap(canvas, yflip=true, axis=nothing, border=:none; kwargs...)
+    _heatmap(canvas, g; kwargs...)
+end
+
+# If the grid has raster dimensions, 
+# plot as a raster on a spatial grid
+function _heatmap(canvas, g; kwargs...)
+    if isnothing(dims(g))
+        heatmap(canvas; yflip=true, axis=nothing, border=:none, aspect_ratio=:equal, kwargs...)
+    else
+        heatmap(Raster(canvas, dims(g)); kwargs...)
+    end
 end
 
 """
diff --git a/src/gridrsp.jl b/src/gridrsp.jl
index 8e9f711..40fc321 100644
--- a/src/gridrsp.jl
+++ b/src/gridrsp.jl
@@ -31,16 +31,27 @@ function GridRSP(g::Grid; θ=nothing)
     return GridRSP(g, θ, Pref, W, Z)
 end
 
+_get_grid(grsp::GridRSP) = grsp.g
+_get_grid(g::Grid)       = g
+
+_maybe_raster(mat::Raster, g) = mat
+_maybe_raster(mat::AbstractMatix, g::Union{Grid,GridRSP}) = 
+    _maybe_raster(mat, dims(g))
+_maybe_raster(mat::AbstractMatix, ::Nothing) = mat
+_maybe_raster(mat::AbstractMatix, dims::Tuple) = Raster(mat, dims)
+
 function Base.show(io::IO, ::MIME"text/plain", grsp::GridRSP)
     print(io, summary(grsp), " of size ", grsp.g.nrows, "x", grsp.g.ncols)
 end
-
 function Base.show(io::IO, ::MIME"text/html", grsp::GridRSP)
     t = string(summary(grsp), " of size ", grsp.g.nrows, "x", grsp.g.ncols)
     write(io, "<h4>$t</h4>")
     show(io, MIME"text/html"(), plot_outdegrees(grsp.g))
 end
 
+DimensionalData.dims(grsp::GridRSP) = dims(grsp.g)
+
+
 """
     betweenness_qweighted(grsp::GridRSP)::Matrix{Float64}
 
@@ -62,7 +73,7 @@ function betweenness_qweighted(grsp::GridRSP)
         bet[grsp.g.id_to_grid_coordinate_list[i]] = v
     end
 
-    return bet
+    return _maybe_raster(bet, grsp)
 end
 
 
@@ -84,7 +95,7 @@ function edge_betweenness_qweighted(grsp::GridRSP)
         [grsp.g.target_qualities[i] for i in grsp.g.id_to_grid_coordinate_list ∩ targetidx],
         targetnodes)
 
-    return betmatrix
+    return _maybe_raster(betmatrix, grsp)
 end
 
 
@@ -139,7 +150,7 @@ function betweenness_kweighted(grsp::GridRSP;
         bet[grsp.g.id_to_grid_coordinate_list[i]] = v
     end
 
-    return bet
+    return _maybe_raster(bet, grsp)
 end
 
 
@@ -359,9 +370,6 @@ function connected_habitat(
 
     return connected_habitat(grsp, S, diagvalue=diagvalue)
 end
-
-_get_grid(grsp::GridRSP) = grsp.g
-_get_grid(g::Grid)       = g
 function connected_habitat(grsp::Union{Grid,GridRSP}, S::Matrix; diagvalue::Union{Nothing,Real}=nothing)
 
     g = _get_grid(grsp)
@@ -383,7 +391,7 @@ function connected_habitat(grsp::Union{Grid,GridRSP}, S::Matrix; diagvalue::Unio
         func[ij] = x
     end
 
-    return func
+    return _maybe_raster(func, grsp)
 end
 
 function connected_habitat(grsp::GridRSP,
@@ -594,5 +602,5 @@ function criticality(grsp::GridRSP;
     landscape = fill(NaN, size(grsp.g))
     landscape[targetidx] = critvec
 
-    return landscape
+    return _maybe_raster(landscape, grsp)
 end
diff --git a/src/utils.jl b/src/utils.jl
index cb6e166..2a9abba 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -9,7 +9,8 @@ function perm_wall_sim(
     corridorpositions=(0.35,0.7),
     impossible_affinity::Real=1e-20,
     nhood_size::Integer=8,
-    kwargs...)
+    kwargs...
+)
 
     # 1. initialize landscape
     affinities = _generate_affinities(nrows, ncols, nhood_size)
@@ -94,7 +95,7 @@ of the cell affinities (costs) weighted by the grid distance (AverageWeight). Th
 respect to eight neighbors (`N8`) or four neighbors (`N4`).
 """
 function graph_matrix_from_raster(
-    R::Matrix;
+    R::AbstractMatrix;
     matrix_type=AffinityMatrix,
     neighbors::Tuple=N8,
     weight=TargetWeight
@@ -173,11 +174,8 @@ function _set_impossible_nodes(g::Grid, node_list::Vector{CartesianIndex{2}}, im
     target_qualities[node_list] .= 0
 
     # Generate a new Grid based on the modified affinities
-    return Grid(size(g)...,
-        affinities=affinities,
-        source_qualities=source_qualities,
-        target_qualities=target_qualities,
-        costs=g.costfunction === nothing ? g.costmatrix : g.costfunction)
+    costs = g.costfunction === nothing ? g.costmatrix : g.costfunction
+    return Grid(size(g)...; affinities, source_qualities, target_qualities, costs)
 end
 
 """

From 498d96f337d1838c068b28625684eb03a81b9621 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Tue, 3 Dec 2024 20:39:43 +0100
Subject: [PATCH 02/11] update examples to use Rasters

---
 examples/1_basic_workflow.jmd            | 90 +++++++++---------------
 examples/2_landmarks.jmd                 | 14 ++--
 examples/3_distance2proximity.jmd        | 69 ++++++++++++------
 examples/4_indep_prob_cost.jmd           | 19 ++---
 examples/5_ConnectedHabitat_variants.jmd | 36 +++++-----
 examples/6_constant_KLd.jmd              | 32 ++++-----
 examples/demo.jmd                        | 10 +--
 src/ConScape.jl                          |  2 +
 src/grid.jl                              | 52 ++++++++------
 src/gridrsp.jl                           |  6 +-
 10 files changed, 168 insertions(+), 162 deletions(-)

diff --git a/examples/1_basic_workflow.jmd b/examples/1_basic_workflow.jmd
index 38bc883..3bcdce4 100644
--- a/examples/1_basic_workflow.jmd
+++ b/examples/1_basic_workflow.jmd
@@ -11,8 +11,9 @@ datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 ```
 
 ```julia
-mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_1000.asc"))
-hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_1000.asc"))
+# TODO remove `replace_missing` and use missingval=NaN in the next Rasters.jl breaking release
+mov_prob = replace_missing(Raster(joinpath(datadir, "mov_prob_1000.asc")), NaN)
+hab_qual = replace_missing(Raster(joinpath(datadir, "hab_qual_1000.asc")), NaN)
 ```
 
 ```julia
@@ -39,12 +40,12 @@ g = ConScape.Grid(size(mov_prob)...,
 ```
 
 ```julia
-ConScape.plot_outdegrees(g, title="Permeabilty")
+plot(ConScape.outdegrees(g); title="Permeabilty")
 # savefig("figure_grid_outdeg.png")
 ```
 
 ```julia
-ConScape.heatmap(g.source_qualities, yflip=true, title="Qualities")
+heatmap(g.source_qualities; title="Qualities")
 # savefig("figure_grid_qualities.png")
 ```
 
@@ -61,7 +62,7 @@ g.affinities
 ```
 
 ```julia
-(g.nrows, g.ncols, g.nrows*g.ncols)
+(g.nrows, g.ncols, g.nrows * g.ncols)
 ```
 
 # Step 2: Habitat creation
@@ -75,7 +76,7 @@ show distance:
 ```julia
 tmp = zeros(5345)
 tmp[4300] = 1
-ConScape.plot_values(g, tmp, title="Target (or is it Source?) Pixel")
+plot(Raster(tmp, g), title="Target (or is it Source?) Pixel")
 ```
 
 ```julia
@@ -83,16 +84,16 @@ dists = ConScape.expected_cost(h);
 ```
 
 ```julia
-ConScape.plot_values(g, dists[:,4300], title="RSP Expected Cost Distance")
+plot(Raster(dists[:,4300], g), title="RSP Expected Cost Distance")
 ```
 
 ```julia
-ConScape.plot_values(g, dists[:,4300], title="RSP Expected Cost Distance")
+plot(Raster(dists[:,4300], g), title="RSP Expected Cost Distance")
 # savefig("figure_ECDistance.png")
 ```
 
 ```julia
-ConScape.plot_values(g, map!(x -> exp(-x/75),dists[:,4300],dists[:,4300]), title="Proximity")
+plot(Raster(map(x -> exp(-x/75), dists[:, 4300]), g), title="Proximity")
 ```
 
 # Step 3: Amount of Connected Habitat
@@ -102,11 +103,11 @@ func = ConScape.connected_habitat(h, distance_transformation=x -> exp(-x/75));
 ```
 
 ```julia
-ConScape.heatmap(Array(func), yflip=true, title="Connected Habitat")
+heatmap(func; title="Connected Habitat")
 ```
 
 ```julia
-ConScape.heatmap(Array(func), yflip=true, title="Connected Habitat")
+heatmap(func; title="Connected Habitat")
 # savefig("figure_ConnectedHabitat.png")
 ```
 
@@ -119,18 +120,18 @@ sum(t -> isfinite(t) ? t : 0.0, func)
 ## quality weighted
 
 ```julia
-ConScape.heatmap(ConScape.betweenness_qweighted(h), yflip=true, title="Quality-weighted Movement Flow")
+plot(ConScape.betweenness_qweighted(h); title="Quality-weighted Movement Flow")
 ```
 
 ```julia
-ConScape.heatmap(ConScape.betweenness_qweighted(h), yflip=true, title="Quality-weighted Movement Flow")
+plot(ConScape.betweenness_qweighted(h); title="Quality-weighted Movement Flow")
 # savefig("figure_Qweighted_flow.png")
 ```
 
 Effect of theta:
 
 ```julia
-tmp = zeros(g.nrows, g.ncols)
+tmp = zeros(dims(g))
 #tmp[42,58] = 1
 #tmp[66, 90] = 1
 tmp[60,70] = 1
@@ -139,64 +140,37 @@ g_tmp = ConScape.Grid(size(mov_prob)...,
                       affinities=ConScape.graph_matrix_from_raster(mov_prob),
                       qualities=tmp,
                       costs=ConScape.MinusLog());
-ConScape.heatmap(g_tmp.source_qualities, yflip=true)
+plot(g_tmp.source_qualities)
 ```
 
 ```julia
 thetas = (2.5, 1.0, 0.5, 0.1, 0.01, 0.001)
-betw = [copy(mov_prob), copy(mov_prob), copy(mov_prob), copy(mov_prob), copy(mov_prob), copy(mov_prob)]
-
-for i in 1:length(thetas)
-    h_tmp = ConScape.GridRSP(g_tmp, θ=thetas[i]);
-    betw[i] = ConScape.betweenness_qweighted(h_tmp)
+betw_rasters = map(thetas) do θ
+    h_tmp = ConScape.GridRSP(g_tmp; θ);
+    ConScape.betweenness_qweighted(h_tmp)
 end
+betw = RasterStack(betw_rasters; name=map(θ -> "theta=$θ", thetas))
 ```
 
 ```julia
-using Plots
-```
-
-```julia
-b1 =ConScape.heatmap(betw[1], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=2.5",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b2 =ConScape.heatmap(betw[2], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=1.0",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b3 =ConScape.heatmap(betw[3], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.5",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b4 =ConScape.heatmap(betw[4], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.1",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b5 =ConScape.heatmap(betw[5], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.01",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b6 =ConScape.heatmap(betw[6], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.001",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-plot(b1,b2,b3,b4,b5,b6, layout = (2,3), size = (2*200, 3*100), dpi=150)
-```
-
-```julia
-b1 =ConScape.heatmap(betw[1], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=2.5",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b2 =ConScape.heatmap(betw[2], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=1.0",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b3 =ConScape.heatmap(betw[3], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.5",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b4 =ConScape.heatmap(betw[4], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.1",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b5 =ConScape.heatmap(betw[5], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.01",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-b6 =ConScape.heatmap(betw[6], yflip=true, legend = :none, xlim=(60,120), ylim=(20,80), title="theta=0.001",
-    titlefont = font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false)
-plot(b1,b2,b3,b4,b5,b6, layout = (2,3), size = (2*200, 3*100), dpi=150)
+plot(betw; xlim=(60, 120), ylim=(20, 80), legend=:none,
+    titlefont=font(10), xaxis=false, yaxis=false, aspect_ratio=:equal, grid=false
+)
 # savefig("output_figures/figure_thetas.png")
 ```
 
 ## Proximity weighted
 
 ```julia
-ConScape.heatmap(ConScape.betweenness_kweighted(h, distance_transformation=x -> exp(-x/75)), yflip=true, title="Proximity-weighted Movement Flow")
+plot(ConScape.betweenness_kweighted(h, distance_transformation=x -> exp(-x/75)); 
+    title="Proximity-weighted Movement Flow"
+)
 ```
 
 ```julia
-ConScape.heatmap(ConScape.betweenness_kweighted(h, distance_transformation=x -> exp(-x/75)), yflip=true, title="Proximity-weighted Movement Flow")
+plot(ConScape.betweenness_kweighted(h, distance_transformation=x -> exp(-x/75)); 
+    title="Proximity-weighted Movement Flow"
+)
 # savefig("output_figures/figure_Pweighted_flow.png")
 ```
 
@@ -207,11 +181,11 @@ ConScape.heatmap(ConScape.betweenness_kweighted(h, distance_transformation=x ->
 We need a smaller (i.e. lower resolution) landscape for computational convenience:
 
 ```julia
-mov_prob_2000, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_2000.asc"))
-hab_qual_2000, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_2000.asc"))
+mov_prob_2000 = replace_missing(Raster(joinpath(datadir, "mov_prob_2000.asc")), NaN)
+hab_qual_2000 = replace_missing(Raster(joinpath(datadir, "hab_qual_2000.asc")), NaN)
 
 # hack for the sensitivity:
-hab_qual_2000[(mov_prob_2000.>0) .== isnan.(hab_qual_2000)] .= 1e-20;
+hab_qual_2000[(mov_prob_2000 .> 0) .== isnan.(hab_qual_2000)] .= 1e-20;
 
 g_2000 = ConScape.Grid(size(mov_prob_2000)...,
                       affinities=ConScape.graph_matrix_from_raster(mov_prob_2000),
diff --git a/examples/2_landmarks.jmd b/examples/2_landmarks.jmd
index 2c8bbdf..975aadd 100644
--- a/examples/2_landmarks.jmd
+++ b/examples/2_landmarks.jmd
@@ -1,7 +1,7 @@
 # Landmark function for testing purposes
 
 ```julia
-using ConScape, SparseArrays, LinearAlgebra
+using ConScape, Rasters, ArchGDAL, SparseArrays, LinearAlgebra, Plots
 ```
 
 ```julia
@@ -19,8 +19,8 @@ datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 ```
 
 ```julia
-mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_1000.asc"))
-hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_1000.asc"));
+mov_prob = replace_missing(Raster(joinpath(datadir, "mov_prob_1000.asc")), NaN)
+hab_qual = replace_missing(Raster(joinpath(datadir, "hab_qual_1000.asc")), NaN);
 ```
 
 ```julia
@@ -78,12 +78,12 @@ cor(filter(!isnan, func), filter(!isnan, func_coarse))
 
 ```julia
 qbetw = ConScape.betweenness_qweighted(h);
-ConScape.heatmap(qbetw, yflip=true)
+heatmap(qbetw)
 ```
 
 ```julia
 qbetw_coarse = ConScape.betweenness_qweighted(h_coarse);
-ConScape.heatmap(qbetw_coarse, yflip=true)
+ConScape.heatmap(qbetw_coarse)
 ```
 
 ```julia
@@ -124,10 +124,6 @@ est_func = first.(tmp)
 cor_func = last.(tmp);
 ```
 
-```julia
-using Plots
-```
-
 ```julia
 sum_func = sum(t -> isnan(t) ? 0.0 : t, func)
 plot(est_func/sum_func)
diff --git a/examples/3_distance2proximity.jmd b/examples/3_distance2proximity.jmd
index f877148..8b0fdc4 100644
--- a/examples/3_distance2proximity.jmd
+++ b/examples/3_distance2proximity.jmd
@@ -1,6 +1,6 @@
 
 ```julia
-using ConScape
+using ConScape, Rasters, ArchGDAL, Plots
 ```
 
 # Step 1: data import and Grid creation
@@ -10,8 +10,8 @@ datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 ```
 
 ```julia
-mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_1000.asc"))
-hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_1000.asc"))
+mov_prob = replace_missing(Raster(joinpath(datadir, "mov_prob_1000.asc")), NaN)
+hab_qual = replace_missing(Raster(joinpath(datadir, "hab_qual_1000.asc")), NaN)
 ```
 
 ```julia
@@ -32,7 +32,7 @@ h = ConScape.GridRSP(g, θ=theta);
 ```julia
 tmp = zeros(5345)
 tmp[4300] = 1
-ConScape.plot_values(g, tmp, title="Target (or is it Source?) Pixel")
+plot(Raster(tmp, g); title="Target (or is it Source?) Pixel")
 ```
 
 # Step 2: Euclidean distances
@@ -42,7 +42,7 @@ euclid = [hypot(xy1[1] - xy2[1], xy1[2] - xy2[2]) for xy1 in g.id_to_grid_coordi
 ```
 
 ```julia
-ConScape.plot_values(g, euclid[:,4300], title="Euclidean Distance")
+plot(Raster(euclid[:, 4300], g); title="Euclidean Distance")
 ```
 
 # Step 3: RSP Expected Cost Distances
@@ -52,49 +52,74 @@ dists = ConScape.expected_cost(h);
 ```
 
 ```julia
-ConScape.plot_values(g, dists[:,4300], title="RSP Expected Cost Distance")
+plot(Raster(dists[:, 4300], g); title="RSP Expected Cost Distance")
 ```
 
 ```julia
-using Plots
-plot(euclid[:,4300], dists[:,4300], seriestype = :scatter, legend=false, xlabel="Euclidean Distance", ylabel="Expected Cost")
+plot(euclid[:, 4300], dists[:, 4300]; 
+    seriestype=:scatter, 
+    legend=false, 
+    xlabel="Euclidean Distance", 
+    ylabel="Expected Cost"
+)
 ```
 
 ```julia
-ConScape.plot_values(g, map!(x -> exp(-x/1000),dists[:,4300],dists[:,4300]), title="RSP Expected Cost Proximity (log)")
+plot(Raster(map(x -> exp(-x / 1000), dists[:,4300]), g); 
+    title="RSP Expected Cost Proximity (log)"
+)
 ```
 
 ```julia
-ConScape.plot_values(g, map!(x -> x<350 ? 1 : 350/x,dists[:,4300],dists[:,4300]), title="RSP Expected Cost Proximity (inv)")
+plot(Raster(map(x -> x<350 ? 1 : 350/x, dists[:,4300]), g); 
+    title="RSP Expected Cost Proximity (inv)"
+)
 ```
 
 ```julia
-plot(map(x -> exp(-x/1000), dists[:,4300]), map(x -> x<350 ? 1 : 350/x, dists[:,4300]), seriestype = :scatter, legend=false, xlabel="log", ylabel="inv")
+plot(
+    map(x -> exp(-x / 1000), dists[:, 4300]), 
+    map(x -> x < 350 ? 1 : 350 / x, dists[:, 4300]); 
+    seriestype=:scatter, 
+    legend=false, 
+    xlabel="log", 
+    ylabel="inv"
+)
 ```
 
 # Step 4: Survival Proximity
 
 ```julia
 surv_prob = h.Z;
-tmp = surv_prob[:,4300]
-map!(t -> t>1 ? 1 : t,  tmp, tmp)
-
-ConScape.plot_values(g, tmp, title="Survival proximity")
+tmp = max.(surv_prob[:, 4300], 1)
+plot(Raster(tmp, g); title="Survival proximity")
 ```
 
 ```julia
-using Plots
 surv_prob = h.Z;
-tmp = surv_prob[:,4300]
-map!(t -> t>1 ? 1 : t,  tmp, tmp)
-
-plot(euclid[:,4300], tmp, seriestype = :scatter, legend=false, xlabel="Euclidean Distance", ylabel="Surival Probability")
+tmp = max.(surv_prob[:, 4300], 1)
+plot(euclid[:, 4300], tmp; 
+    seriestype=:scatter, 
+    legend=false, 
+    xlabel="Euclidean Distance", 
+    ylabel="Surival Probability"
+)
 ```
 
 ```julia
-plot(dists[:,4300], tmp, seriestype = :scatter, legend=false, xlabel="Expected Cost", ylabel="Surival Probability")
+plot(dists[:, 4300], tmp; 
+    seriestype=:scatter, 
+    legend=false, 
+    xlabel="Expected Cost", 
+    ylabel="Surival Probability"
+)
 ```
 
 ```julia
-plot(map(x -> exp(-x/1000), dists[:,4300]), tmp, seriestype = :scatter, legend=false, xlabel="Expected Cost Proximity", ylabel="Surival Probability")
+plot(map(x -> exp(-x / 1000), dists[:, 4300]), tmp; 
+    seriestype=:scatter, 
+    legend=false, 
+    xlabel="Expected Cost Proximity", 
+    ylabel="Surival Probability"
+)
 ```
diff --git a/examples/4_indep_prob_cost.jmd b/examples/4_indep_prob_cost.jmd
index bb8dde1..434f41f 100644
--- a/examples/4_indep_prob_cost.jmd
+++ b/examples/4_indep_prob_cost.jmd
@@ -10,9 +10,9 @@ datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 ```
 
 ```julia
-mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "prob_panther_cropped.asc"))
-mov_cost, meta_c = ConScape.readasc(joinpath(datadir, "mort_panther_cropped.asc"))
-hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "prob_panther_cropped.asc"))
+mov_prob = replace_missing(Raster(joinpath(datadir, "prob_panther_cropped.asc")), NaN)
+mov_cost = replace_missing(Raster(joinpath(datadir, "mort_panther_cropped.asc")), NaN)
+hab_qual = replace_missing(Raster(joinpath(datadir, "prob_panther_cropped.asc")), NaN)
 ```
 
 ```julia
@@ -40,7 +40,7 @@ meta_c
 # Step 2: GridRSP creation
 
 ```julia
-g2 = ConScape.Grid(size(mov_prob)...,
+g2 = ConScape.Grid(size(mov_prob)...;
     affinities=ConScape.graph_matrix_from_raster(mov_cost),
     qualities=hab_qual);
 ConScape.plot_indegrees(g2)
@@ -55,7 +55,7 @@ show survival:
 ```julia
 tmp = zeros(17212)
 tmp[15000] = 1
-ConScape.plot_values(g, tmp, title="Target (or is it Source?) Pixel")
+plot(Raster(tmp, g); title="Target (or is it Source?) Pixel")
 ```
 
 ```julia
@@ -63,13 +63,16 @@ pb = ConScape.survival_probability(h);
 ```
 
 ```julia
-ConScape.plot_values(g, map(t -> t==1 ? NaN : t, pb[:,15000]), title="Survival Probability")
+plot(Raster(map(t -> t==1 ? NaN : t, pb[:,15000]), g); title="Survival Probability")
 ```
 
 ```julia
-func = ConScape.connected_habitat(h, connectivity_function=ConScape.survival_probability, distance_transformation=ConScape.ExpMinus())
+func = ConScape.connected_habitat(h; 
+    connectivity_function=ConScape.survival_probability, 
+    distance_transformation=ConScape.ExpMinus()
+)
 ```
 
 ```julia
-ConScape.heatmap(Array(func), yflip=true, title="Cumulative Connected Habitat (Survival)")
+plot(func; title="Cumulative Connected Habitat (Survival)")
 ```
diff --git a/examples/5_ConnectedHabitat_variants.jmd b/examples/5_ConnectedHabitat_variants.jmd
index e2d9547..8793967 100644
--- a/examples/5_ConnectedHabitat_variants.jmd
+++ b/examples/5_ConnectedHabitat_variants.jmd
@@ -1,6 +1,6 @@
 
 ```julia
-using ConScape
+using ConScape, Rasters, ArchGDAL, Plots
 ```
 
 # Data import and GridRSP creation
@@ -10,8 +10,8 @@ datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 ```
 
 ```julia
-mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_1000.asc"))
-hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_1000.asc"));
+mov_prob = replace_missing(Raster(joinpath(datadir, "mov_prob_1000.asc")), NaN)
+hab_qual = replace_missing(Raster(joinpath(datadir, "hab_qual_1000.asc")), NaN);
 ```
 
 ```julia
@@ -36,11 +36,11 @@ h = ConScape.GridRSP(g, θ=0.001);
 ## Summed Expected Cost
 
 ```julia
-func = ConScape.connected_habitat(h, distance_transformation=x -> exp(-x/2000));
+func = ConScape.connected_habitat(h, distance_transformation=x -> exp(-x / 2000));
 ```
 
 ```julia
-ConScape.heatmap(Array(func), yflip=true)
+plot(func)
 ```
 
 ```julia
@@ -54,11 +54,11 @@ qᵗ = [h.g.target_qualities[i] for i in targetidx]
 ```
 
 ```julia
-similarities = map(t -> iszero(t) ? t : exp(-t/2000), ConScape.expected_cost(h));
+similarities = map(t -> iszero(t) ? t : exp(-t / 2000), ConScape.expected_cost(h));
 ```
 
 ```julia
-ConScape.plot_values(g, similarities[4300,:])
+Raster(similarities[4300, :], g) |> plot
 ```
 
 ```julia
@@ -72,7 +72,7 @@ sum(t -> isnan(t) ? 0.0 : t, func)
 ```
 
 ```julia
-ConScape.plot_values(g, func1)
+Raster(func1, g) |> plot
 ```
 
 ```julia
@@ -87,21 +87,19 @@ landscape = qˢ .* similarities .* qᵗ'
 
 ```julia
 @time vˡ, λ, vʳ = ConScape.eigmax(h, connectivity_function=ConScape.expected_cost, distance_transformation=t -> iszero(t) ? t : exp(-t/2000))
-
 λ
 ```
 
 ```julia
-ConScape.plot_values(g, real.(vʳ))
+Raster(real.(vʳ), g) |> plot
 ```
 
 ```julia
-ConScape.plot_values(g, abs.(real.(vˡ)))
+Raster(abs.(real.(vˡ)), g) |> plot
 ```
 
 ```julia
-using Plots
-plot(func1, real.(vʳ), seriestype = :scatter, legend=false, xlabel="sum", ylabel="eigenvector")
+plot(func1, real.(vʳ); seriestype=:scatter, legend=false, xlabel="sum", ylabel="eigenvector")
 ```
 
 ## Survival
@@ -111,7 +109,7 @@ similarities = h.Z;
 ```
 
 ```julia
-ConScape.plot_values(g, similarities[4300,:])
+Raster(similarities[4300, :], g) |> plot
 ```
 
 ```julia
@@ -121,18 +119,18 @@ func3 = qˢ .* similarities * qᵗ
 ```
 
 ```julia
-ConScape.plot_values(g, func3)
+Raster(func3, g) |> plot
 ```
 
 ```julia
-plot(func1, func3, seriestype = :scatter, legend=false, xlabel="sum", ylabel="eigenvector")
+plot(func1, func3; seriestype=:scatter, legend=false, xlabel="sum", ylabel="eigenvector")
 ```
 
 ## Probability of Connectivity
 
 ```julia
 lcps = ConScape.least_cost_distance(g)
-ConScape.plot_values(g, lcps[4300,:])
+Raster(lcps[4300, :], g) |> plot
 ```
 
 ```julia
@@ -140,7 +138,7 @@ similarities = map(t -> iszero(t) ? t : exp(-t/2.5), lcps);
 ```
 
 ```julia
-ConScape.plot_values(g, similarities[4300,:])
+Raster(similarities[4300, :], g) |> plot
 ```
 
 ```julia
@@ -150,5 +148,5 @@ pc = qˢ .* similarities * qᵗ
 ```
 
 ```julia
-ConScape.plot_values(g, pc)
+Raster(pc, g) |> plot
 ```
diff --git a/examples/6_constant_KLd.jmd b/examples/6_constant_KLd.jmd
index 90144f5..039cb4a 100644
--- a/examples/6_constant_KLd.jmd
+++ b/examples/6_constant_KLd.jmd
@@ -11,8 +11,8 @@ datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 ```
 
 ```julia
-mov_prob, meta_p = ConScape.readasc(joinpath(datadir, "mov_prob_1000.asc"))
-hab_qual, meta_q = ConScape.readasc(joinpath(datadir, "hab_qual_1000.asc"))
+mov_prob = replace_missing(Raster(joinpath(datadir, "mov_prob_1000.asc")), NaN)
+hab_qual = replace_missing(Raster(joinpath(datadir, "hab_qual_1000.asc")), NaN)
 ```
 
 ```julia
@@ -40,7 +40,7 @@ maxkld
 ```
 
 ```julia
-target_mean_D_KL = 0.5*maxkld
+target_mean_D_KL = 0.5 * maxkld
 target_mean_D_KL
 ```
 
@@ -56,7 +56,7 @@ r = optimize(_θ -> mean_D_KL_difference(_θ, g, target_mean_D_KL), 0.0, 3.0; it
 ```
 
 ```julia
-sqrt(r.minimum)/target_mean_D_KL
+sqrt(r.minimum) / target_mean_D_KL
 ```
 
 ```julia
@@ -86,7 +86,7 @@ g = ConScape.perm_wall_sim(nrows, ncols;
     corridorpositions=(0.5,),
     impossible_affinity=0.,
     # Qualities decrease by row
-    qualities=copy(reshape(collect(nrows*ncols:-1:1), ncols, nrows)'))
+    qualities=copy(reshape(collect(nrows * ncols:-1:1), ncols, nrows)'))
 
 
 θ = 10.0
@@ -102,7 +102,7 @@ maxkld
 ```
 
 ```julia
-100*ConScape.mean_kl_divergence(h)/maxkld
+100 * ConScape.mean_kl_divergence(h) / maxkld
 ```
 
 ```julia
@@ -110,9 +110,9 @@ dist = ConScape.expected_cost(h);
 ```
 
 ```julia
-target_node = Int(ceil((ncols-1)/2 - 3)*nrows + ceil((nrows+1)/2))
+target_node = Int(ceil((ncols - 1)/2 - 3) * nrows + ceil((nrows + 1) / 2))
 
-ConScape.plot_values(g, dist[:,target_node])
+plot(Raster(dist[:, target_node], g))
 ```
 
 ```julia
@@ -125,7 +125,7 @@ g2 = ConScape.perm_wall_sim(nrows, ncols;
     corridorpositions=(0.5,),
     impossible_affinity=0.,
     # Qualities decrease by row
-    qualities=copy(reshape(collect(nrows*ncols:-1:1), ncols, nrows)'))
+    qualities=copy(reshape(collect(nrows * ncols:-1:1), ncols, nrows)'))
 ```
 
 ```julia
@@ -149,7 +149,7 @@ r = optimize(_θ -> mean_D_KL_difference(_θ, g2, target_kld), 0.0, 3.0;iteratio
 ```
 
 ```julia
-sqrt(r.minimum)/target_kld
+sqrt(r.minimum) / target_kld
 ```
 
 ```julia
@@ -167,16 +167,16 @@ dist3 = ConScape.expected_cost(h3);
 
 ```julia
 using Plots
-values_orig = ConScape.plot_values(g, map(t -> exp(-t/0.5), dist[:,target_node]))
-Plots.contour(values_orig,fill=true,levels=10, c=:plasma)
+values_orig = ConScape.plot_values(g, map(t -> exp(-t / 0.5), dist[:,target_node]))
+Plots.contour(values_orig; fill=true, levels=10, c=:plasma)
 ```
 
 ```julia
-values_wide = ConScape.plot_values(g2, map(t -> exp(-t/0.5), dist2[:,target_node]))
-Plots.contour(values_wide,fill=true,levels=10, c=:plasma)
+values_wide = ConScape.plot_values(g2, map(t -> exp(-t / 0.5), dist2[:,target_node]))
+Plots.contour(values_wide; fill=true, levels=10, c=:plasma)
 ```
 
 ```julia
-values_consKL = ConScape.plot_values(g2, map(t -> exp(-t/0.5), dist3[:,target_node]))
-Plots.contour(values_consKL,fill=true,levels=10, c=:plasma)
+values_consKL = ConScape.plot_values(g2, map(t -> exp(-t / 0.5), dist3[:,target_node]))
+Plots.contour(values_consKL; fill=true, levels=10, c=:plasma)
 ```
diff --git a/examples/demo.jmd b/examples/demo.jmd
index d6b807e..91e25e8 100644
--- a/examples/demo.jmd
+++ b/examples/demo.jmd
@@ -8,7 +8,7 @@ This is a small demo of the functionalty in the `ConScape` package for Julia.
 First we need to load the package
 
 ```julia
-using ConScape
+using ConScape, Plots
 ```
 
 ```julia
@@ -31,7 +31,9 @@ Cg = ConScape.SimpleWeightedGraph(h.g.costmatrix)
 ```
 
 ```julia
-map!(t -> ifelse(isfinite(t) && !iszero(t), exp(-t), t), distances_all2L_shortestpath, distances_all2L_shortestpath)
+map!(distances_all2L_shortestpath, distances_all2L_shortestpath) do t
+    (isfinite(t) && !iszero(t)) ? exp(-t) : t 
+end
 ```
 
 ```julia
@@ -41,7 +43,7 @@ ConScape.heatmap(reshape(bet_shortest, 30, 60), yflip=true)
 
 ```julia
 bet_q = ConScape.betweenness_qweighted(h)
-ConScape.heatmap(bet_q, yflip=true)
+heatmap(bet_q; yflip=true)
 ```
 
 ```julia
@@ -50,5 +52,5 @@ ConScape.heatmap(bet_q, yflip=true)
 
 ```julia
 bet_k = ConScape.betweenness_kweighted(h);
-ConScape.heatmap(bet_k, yflip=true)
+heatmap(bet_k; yflip=true)
 ```
diff --git a/src/ConScape.jl b/src/ConScape.jl
index 6fb3ea9..3e6b761 100644
--- a/src/ConScape.jl
+++ b/src/ConScape.jl
@@ -2,6 +2,8 @@ module ConScape
 
     using SparseArrays, LinearAlgebra
     using Graphs, Plots, SimpleWeightedGraphs, ProgressLogging, ArnoldiMethod
+    using Rasters
+    using Rasters.DimensionalData
 
     abstract type ConnectivityFunction <: Function end
     abstract type DistanceFunction <: ConnectivityFunction end
diff --git a/src/grid.jl b/src/grid.jl
index f31863e..b1782bd 100644
--- a/src/grid.jl
+++ b/src/grid.jl
@@ -24,9 +24,9 @@ struct Grid
     costfunction::Union{Nothing,Transformation}
     costmatrix::SparseMatrixCSC{Float64,Int}
     id_to_grid_coordinate_list::Vector{CartesianIndex{2}}
-    source_qualities::Matrix{Float64}
+    source_qualities::AbstractMatrix{Float64}
     target_qualities::AbstractMatrix{Float64}
-    dims::Tuple
+    dims::Union{Tuple,Nothing}
 end
 
 """
@@ -63,8 +63,8 @@ function Grid(nrows::Integer,
         throw(ArgumentError("grid size ($nrows, $ncols) is incompatible with size of affinity matrix ($n, $n)"))
     end
 
-    _source_qualities = convert(Matrix{Float64}        , source_qualities)
-    _target_qualities = convert(AbstractMatrix{Float64}, target_qualities)
+    _source_qualities = convert(Matrix{Float64}        , _unwrap(source_qualities))
+    _target_qualities = convert(AbstractMatrix{Float64}, _unwrap(target_qualities))
 
     # Prune
     # id_to_grid_coordinate_list = if prune
@@ -138,7 +138,8 @@ function Base.show(io::IO, ::MIME"text/html", g::Grid)
         write(io, "</td></tr></table>")
     end
 end
-
+_unwrap(R::Raster) = parent(R)
+_unwrap(R::AbstractMatrix) = R
 # Compute a vector of the cartesian indices of nonzero target qualities and
 # the corresponding node id corresponding to the indices
 _targetidx(q::Matrix, grididxs::Vector) = grididxs
@@ -153,32 +154,35 @@ function _targetidx_and_nodes(g::Grid)
     return targetidx, targetnodes
 end
 
-function plot_values(g::Grid, values::Vector; kwargs...)
-    canvas = fill(NaN, g.nrows, g.ncols)
-    for (i,v) in enumerate(values)
-        canvas[g.id_to_grid_coordinate_list[i]] = v
+function _fill_matrix(values, g) 
+    M = fill(NaN, g.nrows, g.ncols)
+    for (i, v) in enumerate(values)
+        M[g.id_to_grid_coordinate_list[i]] = v
     end
-    _heatmap(canvas, g; kwargs...)
+    return M
+end
+
+function Raster(values::AbstractVector, g::Grid; kwargs...)
+    isnothing(dims(g)) && throw(ArgumentError("Grid dims are `nothing` - it was not initialised with a Raster"))
+    return Raster(_fill_matrix(values, g), dims(g); kwargs...)
 end
 
-function plot_outdegrees(g::Grid; kwargs...)
+function outdegrees(g::Grid)
     values = sum(g.affinities, dims=2)
-    canvas = fill(NaN, g.nrows, g.ncols)
-    for (i,v) in enumerate(values)
-        canvas[g.id_to_grid_coordinate_list[i]] = v
-    end
-    _heatmap(canvas, g; kwargs...)
+    _maybe_raster(_fill_matrix(values, g), g)
 end
 
-function plot_indegrees(g::Grid; kwargs...)
+function indegrees(g::Grid; kwargs...)
     values = sum(g.affinities, dims=1)
-    canvas = fill(NaN, g.nrows, g.ncols)
-    for (i,v) in enumerate(values)
-        canvas[g.id_to_grid_coordinate_list[i]] = v
-    end
-    _heatmap(canvas, g; kwargs...)
+    _maybe_raster(_fill_matrix(values, g), g)
 end
 
+plot_values(g::Grid, values::Vector; kwargs...) =
+    _heatmap(_fill_matrix(values, g), g; kwargs...)
+plot_outdegrees(g::Grid; kwargs...) = _heatmap(outdegrees(g), g; kwargs...)
+plot_indegrees(g::Grid; kwargs...) = _heatmap(indegrees(g), g; kwargs...)
+
+
 # If the grid has raster dimensions, 
 # plot as a raster on a spatial grid
 function _heatmap(canvas, g; kwargs...)
@@ -250,7 +254,9 @@ function largest_subgraph(g::Grid)
         g.costfunction === nothing ? g.costmatrix[scci, scci] : mapnz(g.costfunction, affinities),
         g.id_to_grid_coordinate_list[scci],
         g.source_qualities,
-        g.target_qualities)
+        g.target_qualities,
+        g.dims,
+    )
 end
 
 """
diff --git a/src/gridrsp.jl b/src/gridrsp.jl
index 40fc321..351a561 100644
--- a/src/gridrsp.jl
+++ b/src/gridrsp.jl
@@ -35,10 +35,10 @@ _get_grid(grsp::GridRSP) = grsp.g
 _get_grid(g::Grid)       = g
 
 _maybe_raster(mat::Raster, g) = mat
-_maybe_raster(mat::AbstractMatix, g::Union{Grid,GridRSP}) = 
+_maybe_raster(mat::AbstractMatrix, g::Union{Grid,GridRSP}) = 
     _maybe_raster(mat, dims(g))
-_maybe_raster(mat::AbstractMatix, ::Nothing) = mat
-_maybe_raster(mat::AbstractMatix, dims::Tuple) = Raster(mat, dims)
+_maybe_raster(mat::AbstractMatrix, ::Nothing) = mat
+_maybe_raster(mat::AbstractMatrix, dims::Tuple) = Raster(mat, dims)
 
 function Base.show(io::IO, ::MIME"text/plain", grsp::GridRSP)
     print(io, summary(grsp), " of size ", grsp.g.nrows, "x", grsp.g.ncols)

From 423c02980d052b55b52d45b53de45ba0857445b5 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 09:56:33 +0100
Subject: [PATCH 03/11] 1.10 lts as the minimum julia

---
 .github/workflows/CI.yml | 1 -
 Project.toml             | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 0d48154..04fa862 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -16,7 +16,6 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'min'
           - 'lts'
           - '1'
           - 'pre'
diff --git a/Project.toml b/Project.toml
index 274de50..529a8c8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,4 +24,4 @@ Plots = "1.4"
 ProgressLogging = "0.1"
 Rasters = "0.12"
 SimpleWeightedGraphs = "1.1"
-julia = "1.9"
+julia = "1.10"

From f1ab6a2c0c853e5bc51f30da26d36aae5c3e580b Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 09:58:41 +0100
Subject: [PATCH 04/11] somehow CI broke, revert

---
 .github/workflows/CI.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 04fa862..0d48154 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -16,6 +16,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
+          - 'min'
           - 'lts'
           - '1'
           - 'pre'

From 75d4b2e3c22c34a2a59e86fb68d0eb5bfd87f5ad Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 10:07:14 +0100
Subject: [PATCH 05/11] add dev branch to CI

---
 .github/workflows/CI.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 0d48154..92b4544 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,9 +3,11 @@ on:
   pull_request:
     branches:
       - main
+      - dev
   push:
     branches:
       - main
+      - dev
     tags: '*'
   workflow_dispatch:
 jobs:
@@ -16,7 +18,6 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'min'
           - 'lts'
           - '1'
           - 'pre'

From c976c0d147791efe6d2c18ae333dc4e40af0b5d4 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 11:11:59 +0100
Subject: [PATCH 06/11] fix grid construction

---
 src/gridrsp.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/gridrsp.jl b/src/gridrsp.jl
index 351a561..a145bdf 100644
--- a/src/gridrsp.jl
+++ b/src/gridrsp.jl
@@ -431,7 +431,8 @@ function connected_habitat(grsp::GridRSP,
                 grsp.g.costfunction === nothing ? grsp.g.costmatrix : mapnz(grsp.g.costfunction, affinities),
                 grsp.g.id_to_grid_coordinate_list,
                 newsource_qualities,
-                newtarget_qualities)
+                newtarget_qualities,
+                dims(grsp))
 
     newh = GridRSP(newg, θ=grsp.θ)
 

From 76e8d60c577f904d234851b49fd0e1509e3e0d16 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 11:25:31 +0100
Subject: [PATCH 07/11] test Rasters grid

---
 test/runtests.jl | 167 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 165 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index fd32d07..6e69674 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,9 +1,172 @@
 using ConScape, Test, SparseArrays
+using Rasters, ArchGDAL, Plots
 
 datadir = joinpath(dirname(pathof(ConScape)), "..", "data")
 _tempdir = mkdir(tempname())
 
-@testset "sno_2000" begin
+@testset "sno_2000 Rasters" begin
+    landscape = "sno_2000"
+    θ = 0.1
+
+    affinity_raster = Raster(joinpath(datadir, "affinities_$landscape.asc"))
+    affinities = ConScape.graph_matrix_from_raster(affinity_raster)
+    @test affinities[1000:1002, 1000:1002] == [
+        0.0               0.00031508895477488 0.0
+        0.133336775193571 0.0                 0.00119533310704962
+        0.0               0.00031508895477488 0.0]
+
+    qualities = Raster(joinpath(datadir, "qualities_$landscape.asc"))
+    # FIXME! We'd have to handle this somehow in the library
+    @test_broken isnan.(affinity_raster) == isnan.(qualities)
+    qualities[(affinity_raster .> 0) .& isnan.(qualities)] .= 1e-20
+
+    g = ConScape.Grid(size(affinity_raster)...,
+        affinities=affinities,
+        qualities=qualities
+        )
+
+    @testset "Grid fields" begin
+        @test g.ncols == 59
+        @test g.nrows == 44
+        @test g.affinities[1000:1002, 1000:1002] == [
+            0.0               0.557963581550866 0.0
+            0.607917269319296 0.0               0.273319838512689
+            0.0               0.557963581550866 0.0]
+        @test g.id_to_grid_coordinate_list[1000:1002] == [
+            CartesianIndex(37, 41), CartesianIndex(38, 41), CartesianIndex(39, 41)]
+        @test g.source_qualities[30:32, 30:32] == [
+            0.151232712594792 0.146546460358077  0.00748122241586316
+            0.170905506269773 0.0532220626743219 0.00309705330379074
+            0.10284506027383  0.059244283127482  0.0361260830667015]
+        @test g.target_qualities[30:32, 30:32] == [
+            0.151232712594792 0.146546460358077  0.00748122241586316
+            0.170905506269773 0.0532220626743219 0.00309705330379074
+            0.10284506027383  0.059244283127482  0.0361260830667015]
+        @test g.costfunction == ConScape.MinusLog()
+        @test g.costmatrix.nzval[end-2:end] ≈ [
+            17.339919976251554
+            16.99334638597158
+            17.339919976251554]
+        @test dims(g) === dims(affinity_raster)
+    end
+
+    @testset "Rasters are returned" begin
+        @test ConScape.indegrees(g) isa Raster
+        @test ConScape.outdegrees(g)) isa Raster
+        @test Raster(ones(length(g.id_to_grid_coordinate_list)), g) isa Raster
+    end
+
+    grsp = ConScape.GridRSP(g, θ=θ)
+
+    @testset "GridRSP fields" begin
+        @test grsp.Pref.nzval[end-2:end] ≈ [
+            0.00031266870414554466,
+            0.00018055948867712768,
+            0.0001228960184654185]
+        @test grsp.W.nzval[end-2:end] ≈ [
+            5.521044625610329e-5,
+            3.300719810371289e-5,
+            2.1700745653826738e-5]
+        @test grsp.Z[100:102,111:113] ≈ [
+            0.00016367398568399748 0.00015329797258788392 0.00023129137088760695
+            6.432239855543766e-5   5.4944307498337036e-5  7.336349355801132e-5
+            2.421703672476501e-5   1.953276891351377e-5   2.4440692415300965e-5]
+        @test dims(grsp) === dims(affinity_raster)
+    end
+
+    @testset "Test mean_kl_divergence" begin
+        @test ConScape.mean_kl_divergence(grsp) ≈ 323895.3828183995
+    end
+
+    @testset "test adjacency creation with $nn neighbors, $w weighting and $mt" for
+        nn in (ConScape.N4, ConScape.N8),
+            w in (ConScape.TargetWeight, ConScape.AverageWeight),
+                mt in (ConScape.AffinityMatrix, ConScape.CostMatrix)
+# No need to test this on sno_100 and doesn't deepend on θ
+# FIXME! Maybe test mean_kl_divergence for part of the landscape to make sure they all roughly give the same result
+                    @test ConScape.graph_matrix_from_raster(
+                        affinity_raster,
+                        neighbors=nn,
+                        weight=w,
+                        matrix_type=mt) isa ConScape.SparseMatrixCSC
+    end
+
+    @testset "Test betweenness" begin
+        @testset "q-weighted" begin
+            bet = ConScape.betweenness_qweighted(grsp)
+            @test bet isa Raster
+            @test bet[21:23, 21:23] ≈ [
+                1930.1334372152335  256.91061166392745 2866.2998374065373
+                4911.996715311025  1835.991238248377    720.755518530375
+                4641.815380725279  3365.3296878569213   477.1085971945757]
+        end
+
+        @testset "k-weighted" begin
+            bet = ConScape.betweenness_kweighted(grsp, diagvalue=1.)
+            @test bet isa Raster
+            @test bet[21:23, 31:33] ≈ [
+                0.04063917813171917 0.06843246983487516 0.08862506281612659
+                0.03684621201600996 0.10352876485995872 0.1255652231824746
+                0.03190640567704462 0.13832814750469344 0.1961393152256104]
+
+            # Check that summed edge betweennesses corresponds to node betweennesses:
+            bet_edge = ConScape.edge_betweenness_kweighted(grsp, diagvalue=1.)
+            @test bet_edge isa Raster
+            bet_edge_sum = fill(NaN, grsp.g.nrows, grsp.g.ncols)
+            for (i, v) in enumerate(sum(bet_edge,dims=2))
+                bet_edge_sum[grsp.g.id_to_grid_coordinate_list[i]] = v
+            end
+            @test bet_edge_sum[21:23, 31:33] ≈ bet[21:23, 31:33]
+
+            # This is a regression test based on values that we currently believe to be correct
+            bet = ConScape.betweenness_kweighted(grsp, distance_transformation=t -> exp(-t/50))
+            @test bet[21:23, 31:33] ≈ [
+                980.5828087688377 1307.981162399926 1602.8445739784497
+                826.0710054834001 1883.0940077789735 1935.4450344630702
+                676.9212075214159 2228.2700913772774 2884.0409495023364]
+
+            @test ConScape.betweenness_kweighted(grsp, distance_transformation=one)[g.id_to_grid_coordinate_list] ≈
+                ConScape.betweenness_qweighted(grsp)[g.id_to_grid_coordinate_list]
+
+            @test ConScape.edge_betweenness_kweighted(grsp, distance_transformation=one) ≈
+                ConScape.edge_betweenness_qweighted(grsp)
+        end
+    end
+
+    @testset "connected_habitat" begin
+        ch = ConScape.connected_habitat(grsp)
+        @test ch isa Raster{Float64}
+        @test size(ch) == size(grsp.g.source_qualities)
+
+        cl = ConScape.connected_habitat(grsp, CartesianIndex((20,20)))
+        @test cl isa Raster{Float64}
+        @test sum(replace(cl, NaN => 0.0)) ≈ 109.4795495188798
+    end
+
+    @testset "mean_lc_kl_divergence" begin
+        @test ConScape.ConScape.mean_lc_kl_divergence(grsp) ≈ 1.5660600315073947e6
+    end
+
+    @testset "Show methods" begin
+        b = IOBuffer()
+        show(b, "text/plain", g)
+        @test occursin("Grid", String(take!(b)))
+
+        b = IOBuffer()
+        show(b, "text/plain", grsp)
+        @test occursin("GridRSP", String(take!(b)))
+
+        b = IOBuffer()
+        show(b, "text/html", g)
+        @test occursin("Grid", String(take!(b)))
+
+        b = IOBuffer()
+        show(b, "text/html", grsp)
+        @test occursin("GridRSP", String(take!(b)))
+    end
+end
+
+@testset "sno_2000 ascii" begin
     landscape = "sno_2000"
     θ = 0.1
 
@@ -75,7 +238,7 @@ _tempdir = mkdir(tempname())
             17.339919976251554]
     end
 
-    @testset "Grid plotting" begin
+    @testset "Old Grid plotting" begin
         @test ConScape.plot_indegrees(g) isa ConScape.Plots.Plot
         @test ConScape.plot_outdegrees(g) isa ConScape.Plots.Plot
         @test ConScape.plot_values(g,ones(length(g.id_to_grid_coordinate_list))) isa ConScape.Plots.Plot

From a5a187297b208b9701cb2389b707fa9c7d5ca2d6 Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 12:40:49 +0100
Subject: [PATCH 08/11] Add ArchGDAL to test deps

---
 Project.toml | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/Project.toml b/Project.toml
index 529a8c8..0315f9d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,3 +25,9 @@ ProgressLogging = "0.1"
 Rasters = "0.12"
 SimpleWeightedGraphs = "1.1"
 julia = "1.10"
+
+[extras]
+ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+
+[targets]
+test = ["ArchGDAL"]

From 85dca1086b993b8ec56658ff5cea895f8e2e476a Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Thu, 12 Dec 2024 13:34:19 +0100
Subject: [PATCH 09/11] test fixes

---
 test/runtests.jl | 66 ++++++++++++------------------------------------
 1 file changed, 16 insertions(+), 50 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 6e69674..8212ec0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,71 +8,37 @@ _tempdir = mkdir(tempname())
     landscape = "sno_2000"
     θ = 0.1
 
-    affinity_raster = Raster(joinpath(datadir, "affinities_$landscape.asc"))
+    # The way the ascii is read in is reversed and rotated from what GDAL does
+    affinity_raster = reverse(rotr90(replace_missing(Raster(joinpath(datadir, "affinities_$landscape.asc")), NaN)); dims=X)
+    affinity_asc = ConScape.readasc(joinpath(datadir, "affinities_$landscape.asc"))[1]
+    # They are only the same for Float32
+    @test all(Float32.(affinity_asc) .=== Float32.(affinity_raster))
+
     affinities = ConScape.graph_matrix_from_raster(affinity_raster)
-    @test affinities[1000:1002, 1000:1002] == [
+    @test Float32.(affinities[1000:1002, 1000:1002]) == Float32.([
         0.0               0.00031508895477488 0.0
         0.133336775193571 0.0                 0.00119533310704962
-        0.0               0.00031508895477488 0.0]
+        0.0               0.00031508895477488 0.0])
 
-    qualities = Raster(joinpath(datadir, "qualities_$landscape.asc"))
-    # FIXME! We'd have to handle this somehow in the library
-    @test_broken isnan.(affinity_raster) == isnan.(qualities)
-    qualities[(affinity_raster .> 0) .& isnan.(qualities)] .= 1e-20
+    qualities_asc = ConScape.readasc(joinpath(datadir, "qualities_$landscape.asc"))[1]
+    qualities = reverse(rotr90(replace_missing(Raster(joinpath(datadir, "qualities_$landscape.asc")), NaN)); dims=X)
+    @test all(Float32.(qualities_asc) .=== Float32.(qualities))
+    @test dims(qualities) == dims(affinity_raster)
 
     g = ConScape.Grid(size(affinity_raster)...,
         affinities=affinities,
         qualities=qualities
         )
-
-    @testset "Grid fields" begin
-        @test g.ncols == 59
-        @test g.nrows == 44
-        @test g.affinities[1000:1002, 1000:1002] == [
-            0.0               0.557963581550866 0.0
-            0.607917269319296 0.0               0.273319838512689
-            0.0               0.557963581550866 0.0]
-        @test g.id_to_grid_coordinate_list[1000:1002] == [
-            CartesianIndex(37, 41), CartesianIndex(38, 41), CartesianIndex(39, 41)]
-        @test g.source_qualities[30:32, 30:32] == [
-            0.151232712594792 0.146546460358077  0.00748122241586316
-            0.170905506269773 0.0532220626743219 0.00309705330379074
-            0.10284506027383  0.059244283127482  0.0361260830667015]
-        @test g.target_qualities[30:32, 30:32] == [
-            0.151232712594792 0.146546460358077  0.00748122241586316
-            0.170905506269773 0.0532220626743219 0.00309705330379074
-            0.10284506027383  0.059244283127482  0.0361260830667015]
-        @test g.costfunction == ConScape.MinusLog()
-        @test g.costmatrix.nzval[end-2:end] ≈ [
-            17.339919976251554
-            16.99334638597158
-            17.339919976251554]
-        @test dims(g) === dims(affinity_raster)
-    end
+    @test dims(g) === dims(qualities)
 
     @testset "Rasters are returned" begin
         @test ConScape.indegrees(g) isa Raster
-        @test ConScape.outdegrees(g)) isa Raster
+        @test ConScape.outdegrees(g) isa Raster
         @test Raster(ones(length(g.id_to_grid_coordinate_list)), g) isa Raster
     end
 
     grsp = ConScape.GridRSP(g, θ=θ)
-
-    @testset "GridRSP fields" begin
-        @test grsp.Pref.nzval[end-2:end] ≈ [
-            0.00031266870414554466,
-            0.00018055948867712768,
-            0.0001228960184654185]
-        @test grsp.W.nzval[end-2:end] ≈ [
-            5.521044625610329e-5,
-            3.300719810371289e-5,
-            2.1700745653826738e-5]
-        @test grsp.Z[100:102,111:113] ≈ [
-            0.00016367398568399748 0.00015329797258788392 0.00023129137088760695
-            6.432239855543766e-5   5.4944307498337036e-5  7.336349355801132e-5
-            2.421703672476501e-5   1.953276891351377e-5   2.4440692415300965e-5]
-        @test dims(grsp) === dims(affinity_raster)
-    end
+    @test dims(grsp) === dims(affinity_raster)
 
     @testset "Test mean_kl_divergence" begin
         @test ConScape.mean_kl_divergence(grsp) ≈ 323895.3828183995
@@ -95,7 +61,7 @@ _tempdir = mkdir(tempname())
         @testset "q-weighted" begin
             bet = ConScape.betweenness_qweighted(grsp)
             @test bet isa Raster
-            @test bet[21:23, 21:23] ≈ [
+            @test parent(bet[21:23, 21:23]) ≈ [
                 1930.1334372152335  256.91061166392745 2866.2998374065373
                 4911.996715311025  1835.991238248377    720.755518530375
                 4641.815380725279  3365.3296878569213   477.1085971945757]

From 660e36dfe68d849377736621d3fa1c615801d2bd Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 13 Dec 2024 13:24:10 +0100
Subject: [PATCH 10/11] bugfix tests

---
 test/runtests.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 8212ec0..e47a554 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,6 +25,8 @@ _tempdir = mkdir(tempname())
     @test all(Float32.(qualities_asc) .=== Float32.(qualities))
     @test dims(qualities) == dims(affinity_raster)
 
+    qualities[(affinity_raster .> 0) .& isnan.(qualities)] .= 1e-20
+
     g = ConScape.Grid(size(affinity_raster)...,
         affinities=affinities,
         qualities=qualities

From b24287d9cd8663716208bd885d5deecd238a433e Mon Sep 17 00:00:00 2001
From: Rafael Schouten <rafaelschouten@gmail.com>
Date: Fri, 13 Dec 2024 14:23:13 +0100
Subject: [PATCH 11/11] passing tests

---
 src/gridrsp.jl   |  2 +-
 test/runtests.jl | 16 +++++++++-------
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/gridrsp.jl b/src/gridrsp.jl
index a145bdf..6638118 100644
--- a/src/gridrsp.jl
+++ b/src/gridrsp.jl
@@ -95,7 +95,7 @@ function edge_betweenness_qweighted(grsp::GridRSP)
         [grsp.g.target_qualities[i] for i in grsp.g.id_to_grid_coordinate_list ∩ targetidx],
         targetnodes)
 
-    return _maybe_raster(betmatrix, grsp)
+    return betmatrix
 end
 
 
diff --git a/test/runtests.jl b/test/runtests.jl
index e47a554..5abd338 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -63,23 +63,23 @@ _tempdir = mkdir(tempname())
         @testset "q-weighted" begin
             bet = ConScape.betweenness_qweighted(grsp)
             @test bet isa Raster
-            @test parent(bet[21:23, 21:23]) ≈ [
+            @test isapprox(bet[21:23, 21:23], [
                 1930.1334372152335  256.91061166392745 2866.2998374065373
                 4911.996715311025  1835.991238248377    720.755518530375
-                4641.815380725279  3365.3296878569213   477.1085971945757]
+                4641.815380725279  3365.3296878569213   477.1085971945757], atol=1e-3)
         end
 
         @testset "k-weighted" begin
             bet = ConScape.betweenness_kweighted(grsp, diagvalue=1.)
             @test bet isa Raster
-            @test bet[21:23, 31:33] ≈ [
+            @test isapprox(bet[21:23, 31:33], [
                 0.04063917813171917 0.06843246983487516 0.08862506281612659
                 0.03684621201600996 0.10352876485995872 0.1255652231824746
-                0.03190640567704462 0.13832814750469344 0.1961393152256104]
+                0.03190640567704462 0.13832814750469344 0.1961393152256104], atol=1e-6)
 
             # Check that summed edge betweennesses corresponds to node betweennesses:
             bet_edge = ConScape.edge_betweenness_kweighted(grsp, diagvalue=1.)
-            @test bet_edge isa Raster
+            @test bet_edge isa SparseMatrixCSC
             bet_edge_sum = fill(NaN, grsp.g.nrows, grsp.g.ncols)
             for (i, v) in enumerate(sum(bet_edge,dims=2))
                 bet_edge_sum[grsp.g.id_to_grid_coordinate_list[i]] = v
@@ -88,10 +88,12 @@ _tempdir = mkdir(tempname())
 
             # This is a regression test based on values that we currently believe to be correct
             bet = ConScape.betweenness_kweighted(grsp, distance_transformation=t -> exp(-t/50))
-            @test bet[21:23, 31:33] ≈ [
+            # TODO the floating point differnce is more 
+            # significant here, 1e-3 is as gooda as it can get
+            @test isapprox(bet[21:23, 31:33], [
                 980.5828087688377 1307.981162399926 1602.8445739784497
                 826.0710054834001 1883.0940077789735 1935.4450344630702
-                676.9212075214159 2228.2700913772774 2884.0409495023364]
+                676.9212075214159 2228.2700913772774 2884.0409495023364], atol=1e-3)
 
             @test ConScape.betweenness_kweighted(grsp, distance_transformation=one)[g.id_to_grid_coordinate_list] ≈
                 ConScape.betweenness_qweighted(grsp)[g.id_to_grid_coordinate_list]