curiosity-score-bon/01_load_sdms.jl at main · PoisotLab/curiosity-score-bon · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
using SpeciesDistributionToolkit
const SDT = SpeciesDistributionToolkit

using BiodiversityObservationNetworks
const BON = BiodiversityObservationNetworks

using Statistics
import CSV
using CairoMakie

# QC
esri_poly = getpolygon(PolygonData(ESRI, Places))
qc = esri_poly["Place name"=>"Quebec"]

# Get a bbox for QC
bbox = (left = -80., right=-50., bottom=44., top=63.)

# On prend les layers
layers = SDMLayer{Float32}[SDMLayer(raw_sdm; bandnumber=1, bbox...) for raw_sdm in readdir("data/sdm", join=true)]
[clamp!(layer.grid, 0., 1.) for layer in layers]

# Noms des espèces
species = replace.(readdir("data/sdm", join=false), ".tif" => "", "_" => " ")

# Dict pour les SDMs
sdms = Dict(zip(species, layers))

# Curiosity scores
Sc(p; b=1.0) = one(p) - (one(p) - p)^b
Si(p; b=1.0) = ((((p - 0.5)^2.0) - 0.25) / (-0.25))^b

# Visualiser le score
probas = LinRange(0.0, 1.0, 500)
begin
    f = Figure()
    ax = Axis(f[1,1], title="Co-occurrence")
    lines!(ax, probas, Sc.(probas; b=0.1))
    lines!(ax, probas, Sc.(probas; b=1.0))
    lines!(ax, probas, Sc.(probas; b=10.0))
    ax = Axis(f[1,2], title="Interaction")
    lines!(ax, probas, Si.(probas; b=0.1))
    lines!(ax, probas, Si.(probas; b=1.0))
    lines!(ax, probas, Si.(probas; b=10.0))
    f
end

# Load the network
interactions = CSV.File("data/canadian_thresholded_reconciled.csv")
#interactions = filter(x -> x.to == "Sylvilagus floridanus", interactions)

preys = unique([int.to for int in interactions])
preds = unique([int.from for int in interactions])
species_pool = unique(vcat(preys, preds))

# Scores
score_coo = []
score_int = []

# Curiosity score for co-occurrence
for prey in preys
    if prey in keys(sdms)
        for pred in preds
            if pred in keys(sdms)
                pair = sdms[prey] * sdms[pred]
                push!(score_coo, Sc.(pair; b=0.5))
                # Get interaction probability
                int = filter(x -> (x.from == pred) & (x.to == prey), interactions)
                push!(score_int, isempty(int) ? Si(0.0; b=0.1) : Si(only(int).score; b=0.1))
            end
        end
    end
end

# Weighted scores
score = score_coo .* score_int
curiosity = sum(score)


f = Figure()
ax1 = Axis(f[1,1], title="Co-occurrence", aspect=DataAspect())
heatmap!(ax1, sum(score_coo))
ax2 = Axis(f[1,2], title="Curiosity", aspect=DataAspect())
heatmap!(ax2, curiosity)
f

zcurio = (curiosity - mean(curiosity))/std(curiosity)

# Colorbar range
bar_range(x) = maximum(abs.(quantile(x, [0.05, 0.95]))) .* (-1, 1)

# Figure
f = Figure()
ax = Axis(f[1, 1], aspect=DataAspect())
hm = heatmap!(ax, zcurio, colormap=:roma, colorrange=bar_range(zcurio))
lines!(ax, qc, color=:black)
Colorbar(f[1, 2], hm)
hidespines!(ax)
f

# BON

sites = BON.sample(BalancedAcceptance(50), curiosity)

f = Figure()
ax = Axis(f[1, 1], aspect=DataAspect())
hm = heatmap!(ax, curiosity, colormap=:navia)
lines!(ax, qc, color=:black)
Colorbar(f[1, 2], hm)
hidespines!(ax)
scatter!(ax, sites.nodes, color=:transparent, strokewidth=2, strokecolor=:white)
f