Merge branch 'feat/ssl-update' into dev

Project.toml

@@ -18,3 +18,6 @@ SimpleSDMLayers = "2c645270-77db-11e9-22c3-0f302a89c64c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
+
+[compat]
+julia = "1.5"

data/proc/glossary.csv

@@ -34,10 +34,10 @@ sp3,species,Setophaga_petechia,Yellow Warbler
 sp4,species,Setophaga_ruticilla,American Redstart
 sp5,species,Mniotilta_varia,Black-and-white Warbler
 sp6,species,Oreothlypis_celata,Orange-crowned Warbler
-sp7,species,Seiurus_aurocapilla,Ovenbird
-sp8,species,Setophaga_palmarum,Palm Warbler
-sp9,species,Setophaga_pinus,Pine Warbler
-sp10,species,Setophaga_americana,Northern Parula
+sp7,species,Setophaga_palmarum,Palm Warbler
+sp8,species,Seiurus_aurocapilla,Ovenbird
+sp9,species,Setophaga_americana,Northern Parula
+sp10,species,Setophaga_pinus,Pine Warbler
 sp11,species,Setophaga_virens,Black-throated Green Warbler
 sp12,species,Cardellina_pusilla,Wilson's Warbler
 sp13,species,Setophaga_magnolia,Magnolia Warbler
@@ -53,12 +53,12 @@ sp22,species,Setophaga_discolor,Prairie Warbler
 sp23,species,Setophaga_citrina,Hooded Warbler
 sp24,species,Setophaga_townsendi,Townsend's Warbler
 sp25,species,Vermivora_cyanoptera,Blue-winged Warbler
-sp26,species,Parkesia_motacilla,Louisiana Waterthrush
-sp27,species,Protonotaria_citrea,Prothonotary Warbler
-sp28,species,Setophaga_tigrina,Cape May Warbler
+sp26,species,Setophaga_tigrina,Cape May Warbler
+sp27,species,Parkesia_motacilla,Louisiana Waterthrush
+sp28,species,Protonotaria_citrea,Prothonotary Warbler
 sp29,species,Setophaga_nigrescens,Black-throated Gray Warbler
-sp30,species,Cardellina_canadensis,Canada Warbler
-sp31,species,Setophaga_castanea,Bay-breasted Warbler
+sp30,species,Setophaga_castanea,Bay-breasted Warbler
+sp31,species,Cardellina_canadensis,Canada Warbler
 sp32,species,Geothlypis_tolmiei,MacGillivray's Warbler
 sp33,species,Helmitheros_vermivorum,Worm-eating Warbler
 sp34,species,Geothlypis_philadelphia,Mourning Warbler
@@ -71,23 +71,22 @@ sp40,species,Myioborus_pictus,Painted Redstart
 sp41,species,Oreothlypis_virginiae,Virginia's Warbler
 sp42,species,Setophaga_graciae,Grace's Warbler
 sp43,species,Limnothlypis_swainsonii,Swainson's Warbler
-sp44,species,Oporornis_agilis,Connecticut Warbler
-sp45,species,Basileuterus_rufifrons,Rufous-capped Warbler
-sp46,species,Cardellina_rubrifrons,Red-faced Warbler
-sp47,species,Setophaga_chrysoparia,Golden-cheeked Warbler
-sp48,species,Myioborus_miniatus,Slate-throated Redstart
-sp49,species,Setophaga_pitiayumi,Tropical Parula
-sp50,species,Basileuterus_culicivorus,Golden-crowned Warbler
-sp51,species,Oreothlypis_superciliosa,Crescent-chested Warbler
-sp52,species,Geothlypis_poliocephala,Gray-crowned Yellowthroat
-sp53,species,Setophaga_kirtlandii,Kirtland's Warbler
-sp54,species,Basileuterus_belli,Golden-browed Warbler
-sp55,species,Cardellina_rubra,Red Warbler
-sp56,species,Basileuterus_lachrymosus,Fan-tailed Warbler
-sp57,species,Oreothlypis_crissalis,Colima Warbler
-sp58,species,Geothlypis_beldingi,Belding's Yellowthroat
-sp59,species,Geothlypis_nelsoni,Hooded Yellowthroat
-sp60,species,Cardellina_versicolor,Pink-headed Warbler
+sp44,species,Cardellina_rubrifrons,Red-faced Warbler
+sp45,species,Oporornis_agilis,Connecticut Warbler
+sp46,species,Setophaga_chrysoparia,Golden-cheeked Warbler
+sp47,species,Setophaga_pitiayumi,Tropical Parula
+sp48,species,Basileuterus_rufifrons,Rufous-capped Warbler
+sp49,species,Setophaga_kirtlandii,Kirtland's Warbler
+sp50,species,Myioborus_miniatus,Slate-throated Redstart
+sp51,species,Geothlypis_poliocephala,Gray-crowned Yellowthroat
+sp52,species,Basileuterus_culicivorus,Golden-crowned Warbler
+sp53,species,Oreothlypis_superciliosa,Crescent-chested Warbler
+sp54,species,Basileuterus_lachrymosus,Fan-tailed Warbler
+sp55,species,Oreothlypis_crissalis,Colima Warbler
+sp56,species,Geothlypis_beldingi,Belding's Yellowthroat
+sp57,species,Basileuterus_belli,Golden-browed Warbler
+sp58,species,Cardellina_rubra,Red Warbler
+sp59,species,Geothlypis_flavovelata,Altamira Yellowthroat
+sp60,species,Geothlypis_nelsoni,Hooded Yellowthroat
 sp61,species,Geothlypis_speciosa,Black-polled Yellowthroat
-sp62,species,Geothlypis_flavovelata,Altamira Yellowthroat
-sp63,species,Vermivora_bachmanii,Bachman's Warbler
+sp62,species,Vermivora_bachmanii,Bachman's Warbler

src/00b_preparation_ebd-preparation.jl

@@ -9,7 +9,7 @@ end
 ## EBD data preparation
 
 # Load data from CSV files (from file cut with terminal)
-@time df = DataFrame!(CSV.File(joinpath("data", "raw", "ebd_warblers_cut.csv"), header=true, delim="\t", copycols = true))
+@time df = DataFrame!(CSV.File(joinpath("data", "raw", "ebd_warblers_cut.csv"), header=true, delim="\t"))
 
 # Prepare data (arrange values & columns)
 @time prepare_ebd_data!(df)

src/00c_preparation_landcover.jl

@@ -34,18 +34,18 @@ end
 coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
 
 # Test loading variables
-lc_vars = landcover(1:10, resolution = "5")
-lc_vars = landcover(1:10, resolution = "10")
-lc_vars = map(x -> landcover(x, resolution = "10")[coords], 1:10)
+lc_vars = landcover(1:10, resolution = 5.0)
+lc_vars = landcover(1:10, resolution = 10.0)
+lc_vars = map(x -> landcover(x, resolution = 10.0)[coords], 1:10)
 
 # Load worldclim variables to compare
-wc_vars = map(x -> worldclim(x, resolution = "10")[coords], 1:19);
+wc_vars = map(x -> worldclim(x, resolution = 10.0)[coords], 1:19);
 
 ## Plot environmental variables examples
 # Plot wcvars1 (temperature)
 wc_plot = plotSDM2(wc_vars[1], colorbar_title = "Annual Mean Temperature (°C)")
 # Plot lcvars2 (urban)
-lc_plot = plotSDM(lc_vars[2], colorbar_title = "Crops land cover (%)")
+lc_plot = plotSDM2(lc_vars[2], colorbar_title = "Crops land cover (%)")
 
 # Get variable info
 glossary = DataFrame!(CSV.File(joinpath("data", "proc", "glossary.csv")))
@@ -93,6 +93,7 @@ end
 env_vars = [wc_vars; lc_vars]
 # Create env matrix
 env_mat = mapreduce(x -> vec(x.grid), hcat, env_vars)
+replace!(x -> isnothing(x) ? NaN : x, env_mat)
 # Create env dataframe
 env_df = DataFrame(env_mat)
 rename!(env_df, vcat(Symbol.("wc", 1:size(wc_vars, 1)), Symbol.("lc", 1:size(lc_vars, 1))))

src/01_distributions.jl

@@ -21,32 +21,46 @@ else
 end
 
 ## Get & prepare data
+# Define coordinates range
+coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
+copy_layer = worldclim(1)[coords]
+coords_bbox = (left = copy_layer.left, right = copy_layer.right, 
+               bottom = copy_layer.bottom, top = coords.top)
 # Load data from CSV files
-@time df = DataFrame!(CSV.File(joinpath("data", "proc", "ebd_warblers_prep.csv"), header=true, delim="\t"))
+df = DataFrame!(CSV.File(joinpath("data", "proc", "ebd_warblers_prep.csv"), header=true, delim="\t"))
+# Remove groupIdentifier column (causing bug)
+select!(df, Not(:groupIdentifier))
+# Filter observations outside coordinates range
+filter!(x -> coords_bbox.left < x.longitude < coords_bbox.right, df)
+filter!(x -> coords_bbox.bottom < x.latitude < coords_bbox.top, df)
 # Separate species
-warblers = [df[df.species .== u,:] for u in unique(df.species)]
+# warblers = [df[df.species .== u,:] for u in unique(df.species)]
+warblers = groupby(df, :species)
 # Reorder species by frequency
-sort!(warblers, by = x -> nrow(x), rev = true)
+# sort!(warblers, by = x -> nrow(x), rev = true)
+warblers = warblers |> 
+    x -> combine(nrow, x) |>
+    x -> sortperm(x, :nrow, rev = true) |>
+    x -> warblers[x]
+warblers = [w for w in warblers]
 # Extract species names
 spenames = [w.species[1] for w in warblers]
 specommon = [w.commonName[1] for w in warblers]
 # Create index Dict for species names
 speindex = indexmap(spenames)
 
-# Define coordinates range
-coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
 # Observed coordinates range
 coords_obs = (left = minimum(df.longitude), right = maximum(df.longitude),
               bottom = minimum(df.latitude), top = maximum(df.latitude))
 
 ## Get environmental data (with different training resolutions)
 # WorldClim data
-wc_vars = map(x -> worldclim(x, resolution = "10")[coords], [1,12]);
+wc_vars = map(x -> worldclim(x, resolution = 10.0)[coords], [1,12]);
 # Landcover data
-lc_vars = map(x -> landcover(x, resolution = "10")[coords], 1:10);
+lc_vars = map(x -> landcover(x, resolution = 10.0)[coords], 1:10);
 # Training data with finer resolution
-wc_vars_train = map(x -> worldclim(x, resolution = "5")[coords_obs], [1,12]);
-lc_vars_train = map(x -> landcover(x, resolution = "5")[coords_obs], 1:10);
+wc_vars_train = map(x -> worldclim(x, resolution = 5.0)[coords_obs], [1,12]);
+lc_vars_train = map(x -> landcover(x, resolution = 5.0)[coords_obs], 1:10);
 
 # Combine environmental data
 env_vars = vcat(wc_vars, lc_vars)
@@ -60,7 +74,7 @@ if (@isdefined create_distributions) && create_distributions == true
     # Select function to run given desired outcome
     if outcome == "raw"
         # Get raw distributions
-        @time distributions = @showprogress pmap(x -> presence_absence(x, env_vars[1]), warblers)
+        @time distributions = @showprogress map(x -> presence_absence(x, env_vars[1]), warblers)
         # @time distributions = @showprogress pmap(x -> presence_absence(x, env_vars_train[1], full_range = true, binary = false), warblers)
     elseif outcome == "bio"
         # Get sdm distributions (with different training resolutions)
@@ -101,24 +115,24 @@ else
 end
 
 ## Count sites with presence per species
-pres_counts = [length(filter(x -> x > 0.0, species.grid)) for species in distributions]
+pres_counts = [length(filter(x -> !isnothing(x) && x > 0.0, species.grid)) for species in distributions]
 sort(pres_counts)
 
 ## Plot result
 # Species 1
 sp1 = "Setophaga_townsendi"
-map_sp1 = plotSDM(distributions[speindex[sp1]], c=:BuPu,
-                  title = "$(sp1) distribution ($(outcome))",
-                  colorbar=:none, dpi=300)
-scatter!(map_sp1, [NaN], label="Occurrence", color=:purple, markershape=:rect, markersize=2,
-                        legend=:bottomright, legendfontsize=5)
+map_sp1 = plotSDM2(distributions[speindex[sp1]], c = reverse(cgrad(:BuPu)),
+                   title = "$(sp1) distribution ($(outcome))",
+                   colorbar = :none)
+scatter!(map_sp1, [NaN], label = "Occurrence", color = :purple, markershape = :rect, markersize = 2,
+                        legend = :bottomright, legendfontsize = 5)
 # Species 2
 sp2 = "Setophaga_petechia"
-map_sp2 = plotSDM(distributions[speindex[sp2]], c=:BuPu,
-                  title = "$(sp2) distribution ($(outcome))",
-                  colorbar=:none, dpi=300)
-scatter!(map_sp2, [NaN], label="Occurrence", color=:purple, markershape=:rect, markersize=2,
-                        legend=:bottomright, legendfontsize=5)
+map_sp2 = plotSDM2(distributions[speindex[sp2]], c = reverse(cgrad(:BuPu)),
+                   title = "$(sp2) distribution ($(outcome))",
+                   colorbar = :none)
+scatter!(map_sp2, [NaN], label = "Occurrence", color = :purple, markershape = :rect, markersize = 2,
+                        legend = :bottomright, legendfontsize = 5)
 
 ## Export figures
 # save_figures = true # should figures be overwritten (optional)

src/02c_training_bart.R

@@ -87,8 +87,13 @@ for (gp in seq_along(spe_groups)) {
 
     message(paste0("Training multi-species group (", gp, "/", length(spe_groups)), ")")
 
+    ## Create BART models
     # create_models <- TRUE
+    modelname <- ifelse(exists("subset_qc") && isTrue(subset_qc), paste0("bart_models_qc", gp, ".RData"), paste0("bart_models", gp, ".RData"))
+    filepath <- here("data", "rdata", modelname)
     if (exists("create_models") && isTRUE(create_models)){
+        # Run models
+        message("Creating models in parallel: ", modelname)
         set.seed(42)
         system.time(
             sdms <- future_map(
@@ -100,16 +105,16 @@ for (gp in seq_along(spe_groups)) {
                 .progress = TRUE
             )
         ) # ~ 4 min in parallel
-    }
-
-    # Export results
-    # save_models <- TRUE
-    filepath <- here("data", "rdata", paste0("bart_models", gp, ".RData"))
-    if (exists("save_models") && isTRUE(save_models)) {
-        message("Saving models to RData")
-        save(sdms, file = filepath)
+
+        # Export results
+        # save_models <- TRUE
+        if (exists("save_models") && isTRUE(save_models)) {
+            message("Saving models to RData: ", modelname)
+            save(sdms, file = filepath)
+        }
     } else {
-        message("Loading models from RData")
+        # Load models from files
+        message("Loading models from RData: ", modelname)
         load(filepath)
     }
 
@@ -143,7 +148,7 @@ for (gp in seq_along(spe_groups)) {
 
     ## 4. Multi-species predictions ####
 
-    message("Predicting species distributions")
+    message("Predicting species distributions: ", modelname)
 
     # Quantile Predictions
     system.time(
@@ -208,6 +213,14 @@ for (gp in seq_along(spe_groups)) {
 }
 )
 
+# Check results
+results_global
+varimps_global
+pred_df_global
+lower_df_global
+upper_df_global
+pres_df_global
+
 # Export to CSV
 # save_predictions <- TRUE
 if (exists("save_predictions") && isTRUE(save_predictions)) {
@@ -244,25 +257,28 @@ spe_sel <- c("sp1", "sp6", "sp17")
 vars_sel <- map(spe_sel, ~ NULL)
 names(vars_sel) <- spe_sel
 # Run variable selection
-tictoc::tic("total")
-for(sp in spe_sel){
-    tictoc::tic(sp)
-    set.seed(42)
-    message(paste0("Variable selection for ", sp, " (", which(sp == spe_sel), "/", length(spe_sel)), ")")
-    # Save plot to png
-    png(here("fig", "bart", paste0("x_bart_vars-select_", sp, ".png")))
-    step_vars <- variable.step(
-        y.data = spe[[sp]], 
-        x.data = env[xnames],
-        iter = 50
-    )
-    dev.off()
-    # Save variables to list
-    vars_sel[[sp]] <- step_vars
+# variable_selection <- TRUE
+if (exists("variable_selection") && isTRUE(variable_selection)) {
+    tictoc::tic("total")
+    for(sp in spe_sel){
+        tictoc::tic(sp)
+        set.seed(42)
+        message(paste0("Variable selection for ", sp, " (", which(sp == spe_sel), "/", length(spe_sel)), ")")
+        # Save plot to png
+        png(here("fig", "bart", paste0("x_bart_vars-select_", sp, ".png")))
+        step_vars <- variable.step(
+            y.data = spe[[sp]], 
+            x.data = env[xnames],
+            iter = 50
+        )
+        dev.off()
+        # Save variables to list
+        vars_sel[[sp]] <- step_vars
+        tictoc::toc()
+    }
     tictoc::toc()
+    vars_sel
 }
-tictoc::toc()
-vars_sel
 
 ## 3 species QC scale
 # $sp1

src/03b_predictions_random-forests.jl

@@ -56,10 +56,10 @@ end
 ## Create Y matrices
 
 # Get matrix Y
-Y = replace(predictions_full, missing => NaN)
-# Set values to NaN if no species present
-inds_zeros = findall(map(x -> all(iszero.(x)), eachrow(Y)))
-Y[inds_zeros,:] .= NaN
+Y = replace(predictions_full, missing => nothing) |> Array{Union{Nothing, Float32}}
+# Set values to nothing if no species present
+inds_zeros = _indsnotobs(Y)
+Y[inds_zeros, :] .= nothing
 
 ## Create distributions
 
@@ -72,9 +72,9 @@ lims = (left = raw_distributions[1].left, right = raw_distributions[1].right,
         bottom = raw_distributions[1].bottom, top = raw_distributions[1].top)
 
 # Create RF distribution layers
-Ydistrib = replace(Y, 0.0 => NaN)
+Ydistrib = replace(Y, 0.0 => nothing)
 Ygrids = [Ydistrib[:, col] for col in 1:size(Ydistrib,2)]
-Ygrids = reshape.(Ygrids, dims...)
+Ygrids = reshape.(Ygrids, dims...) .|> Array
 distributions = SimpleSDMResponse.(Ygrids, lims...)
 
 ## Export results

src/03c_predictions_bart.jl

@@ -20,13 +20,13 @@ pres_df  = DataFrame!(CSV.File(joinpath("data", "proc", "bart_predictions_pres.c
 ## Create Y matrices
 
 # Get matrix Y
-Y = replace(Array(pres_df), missing => NaN)
-Yprob  = replace(Array(pred_df),  missing => NaN)
-Ylower = replace(Array(lower_df), missing => NaN)
-Yupper = replace(Array(upper_df), missing => NaN)
-# Set values to NaN if no species present
-inds_zeros = findall(map(x -> all(iszero.(x)), eachrow(Y)))
-Y[inds_zeros,:] .= NaN
+Y = replace(Array(pres_df), missing => nothing) |> Array{Union{Nothing, Float32}}
+Yprob  = replace(Array(pred_df),  missing => nothing) |> Array{Union{Nothing, Float32}}
+Ylower = replace(Array(lower_df), missing => nothing) |> Array{Union{Nothing, Float32}}
+Yupper = replace(Array(upper_df), missing => nothing) |> Array{Union{Nothing, Float32}}
+# Set values to nothing if no species present
+inds_zeros = _indsnotobs(Y)
+Y[inds_zeros,:] .= nothing
 
 ## Create distributions
 
@@ -44,9 +44,9 @@ lims = (left = raw_distributions[1].left, right = raw_distributions[1].right,
 # Create distribution layers
 layers = []
 for Y in (Y, Yprob, Ylower, Yupper)
-    Ydistrib = replace(Y, 0.0 => NaN)
+    Ydistrib = replace(Y, 0.0 => nothing)
     Ygrids = [Ydistrib[:, col] for col in 1:size(Ydistrib,2)]
-    Ygrids = reshape.(Ygrids, dims...)
+    Ygrids = reshape.(Ygrids, dims...) .|> Array
     distributions = SimpleSDMResponse.(Ygrids, lims...)
     push!(layers, distributions)
 end