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' diff --git a/Project.toml b/Project.toml index 872d0e3..0315f9d 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,12 @@ Graphs = "1" LaTeXStrings = "1.1" Plots = "1.4" ProgressLogging = "0.1" +Rasters = "0.12" SimpleWeightedGraphs = "1.1" -julia = "1.9" +julia = "1.10" + +[extras] +ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" + +[targets] +test = ["ArchGDAL"] 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 e79e676..b1782bd 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -24,8 +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::Union{Tuple,Nothing} 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) @@ -62,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 @@ -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) @@ -135,7 +138,8 @@ function Base.show(io::IO, ::MIME"text/html", g::Grid) write(io, "") 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 @@ -150,30 +154,43 @@ 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, yflip=true, axis=nothing, border=:none, aspect_ratio=:equal; kwargs...) + return M end -function plot_outdegrees(g::Grid; kwargs...) +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 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, yflip=true, axis=nothing, border=:none; 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 + _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...) + if isnothing(dims(g)) + heatmap(canvas; yflip=true, axis=nothing, border=:none, aspect_ratio=:equal, kwargs...) + else + heatmap(Raster(canvas, dims(g)); kwargs...) end - heatmap(canvas, yflip=true, axis=nothing, border=:none; kwargs...) end """ @@ -237,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 8e9f711..6638118 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::AbstractMatrix, g::Union{Grid,GridRSP}) = + _maybe_raster(mat, dims(g)) +_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) 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, "

$t

") 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 @@ -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, @@ -423,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.θ) @@ -594,5 +603,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 """ diff --git a/test/runtests.jl b/test/runtests.jl index fd32d07..5abd338 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,142 @@ 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 + + # 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 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]) + + 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) + + qualities[(affinity_raster .> 0) .& isnan.(qualities)] .= 1e-20 + + g = ConScape.Grid(size(affinity_raster)..., + affinities=affinities, + qualities=qualities + ) + @test dims(g) === dims(qualities) + + @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, θ=θ) + @test dims(grsp) === dims(affinity_raster) + + @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 isapprox(bet[21:23, 21:23], [ + 1930.1334372152335 256.91061166392745 2866.2998374065373 + 4911.996715311025 1835.991238248377 720.755518530375 + 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 isapprox(bet[21:23, 31:33], [ + 0.04063917813171917 0.06843246983487516 0.08862506281612659 + 0.03684621201600996 0.10352876485995872 0.1255652231824746 + 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 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 + 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)) + # 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], 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] + + @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 +208,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