Setup
Load packages and functions for processing data. The koala_data_transect_view() and koala_data_presence_view() takes the data from proofsafe, as well as the transect spatial data and combines them to cleaned and curated datasets.
Show code
write_data <- FALSE
write_evcs <- FALSE
ala_write <- FALSE
gg_nocturnal_write <- FALSE
koala_translocations_write <- FALSE
knitr::opts_chunk$set(echo = TRUE,
warning = FALSE,
message = FALSE)
mapview::mapviewOptions(fgb = FALSE)
library(sf)
library(dplyr)
library(lidR)
library(rgl)
library(terra)
library(lidRviewer)
library(mapview)
library(galah)
library(DBI)
library(VicmapR)
library(kableExtra)
library(data.table)
library(terra)
library(tidyterra)
library(units)
library(ggplot2)
library(dismo)
library(SDMtune)
library(purrr)
library(cmdstanr)
library(bayesplot)
library(ggspatial)
library(viridis)
library(gridExtra)
library(MASS)
options(mc.cores=4)
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine(override = FALSE)
galah_config(email = "justin.cally@delwp.vic.gov.au", caching = TRUE)
delwp_cols <- c(Teal = "#00B2A9",
Navy = "#201547",
Environment = "#CEDC00",
`Climate Change` = "#FDDA24",
Energy = "#0072CE",
Water = "#71C5E8",
Planning = "#642667",
Infrastructure = "#AF272F",
FFR = "#E57200",
Corporate = "#201547")
delwp_cols2 <- c(delwp_cols, "#737373", "#00B2A9")
gips_bound <- vicmap_query("datavic:VMLITE_LGA") %>%
filter(LGA_NAME == "EAST GIPPSLAND") %>%
collect()
vic_bound <- vicmap_query("datavic:VMLITE_VICTORIA_POLYGON_SU5") %>%
filter(STATE == "VIC") %>% # Philip Island and French Island
collect() %>%
st_transform(3111)
evc_sub <- vicmap_query("datavic:FLORAFAUNA1_VBIOREG100") %>%
filter(INTERSECTS(gips_bound)) %>%
collect() %>%
sf::st_intersection(gips_bound)
delwp_region <- vicmap_query("datavic:VMADMIN_DELWP_REGION") %>%
collect()%>%
st_transform(3111) %>%
sf::st_simplify(dTolerance = 100)
con <- RPostgreSQL::dbConnect(RPostgres::Postgres(),
host = '10.110.7.201',
dbname = 'ari-dev-weda-psql-01',
user = "psql_admin",
password = keyring::key_get(service = "ari-dev-weda-psql-01", username = "psql_admin"),
port = 5432,
service = NULL,
list(sslmode = "require"))
#### Load Existing Survey Data ####
compiled_single_counts <- sf::st_read(con, layer = Id(schema = "koala", table = "compiled_single_counts"))
compiled_double_counts <- sf::st_read(con, layer = Id(schema = "koala", table = "compiled_double_counts"))
#### FUNCTIONS ####
# Transect data
koala_data_transect_view <- function(raw_headers, koala_transects) {
transect_info <- raw_headers %>%
dplyr::mutate(Observer = dplyr::case_when(H0_Observer == "Other" ~ H0_ObserverOther,
TRUE ~ H0_Observer)) %>%
dplyr::mutate(Method = dplyr::case_when(Observer == "Thermal" | stringr::str_detect(Observer, "Thermal") ~ "Thermal",
Observer == "Spotlight" | stringr::str_detect(Observer, "Spot") ~ "Spotlight",
TRUE ~ "Double Observer"),
EffortTime = H1_End_time-H0_Start_time) %>%
dplyr::select(Data_File_Id,
User_Name,
Site = H0_SiteID,
Transect = H0_Transect,
Method,
StartDate = H0_Date,
StartTime = H0_Start_time,
Observer,
ObserverPosition = H0_Observer1,
Temperature = H0_Temp_C,
EndTime = H1_End_time,
EffortTime,
Notes = H1_Notes) %>%
dplyr::left_join(koala_transects %>%
dplyr::mutate(Site = as.numeric(Site),
Transect = as.numeric(Transect)),
by = c("Site", "Transect"))
transect_info
}
# Animal data
koala_data_presence_view <- function(raw_headers, raw_animals, koala_transects) {
raw_animals_simp <- raw_animals %>%
dplyr::mutate(Species = dplyr::case_when(I0_Species == "Other" ~ I0_Animal_sp_other,
TRUE ~ I0_Species),
TreeSpecies = dplyr::case_when(I0_TreeSpecies == "Other" ~ as.character(I0_Tree_sp_other),
TRUE ~ I0_TreeSpecies),
Height = sqrt(I0_DistToAnimal^2 - I0_DistToTree^2)) %>%
dplyr::select(Data_File_Id,
Animal_number = I0_Animal,
Species,
ObservationTime = I0_AnimObsTime,
ObservationType = I0_SeenHeard,
Direction = I0_LorRtrans,
Waypoint = I0_Waypoint_no_,
Latitude = I0_Latitude,
Longitude = I0_Longitude,
HorizontalDistance = I0_DistToTree,
Distance = I0_DistToAnimal,
Angle = I0_AngleToAnimal,
Height,
Bearing = I0_BearingToAnimal,
DistanceAlongTransect = I0_Dist_F_Transect,
TreeSpecies,
SeenByOther = I0_SeenX2,
VBAUpload = I0_VBAUse,
Comments = I0_Comments
)
animals_trans <- koala_data_transect_view(raw_headers, koala_transects) %>%
right_join(raw_animals_simp, by = "Data_File_Id")
animals_trans
}
# Extract EVCs
extract_evc.vicmap_promise <- function(x = VicmapR::vicmap_query("datavic:FLORAFAUNA1_NV2005_EVCBCS"),
y,
group_vars = c("BIOREGION", "X_GROUPNAME"),
EVC_local = NULL,
...) {
# chunk limit
o <- getOption("vicmap.chunk_limit")
options(vicmap.chunk_limit = 1500)
if(is.null(EVC_local)) {
# Pull data from geoserver
EVCs <- x %>%
VicmapR::select(!!!group_vars) %>%
VicmapR::filter(VicmapR::INTERSECTS(y)) %>%
VicmapR::collect() %>%
sf::st_transform(sf::st_crs(y))
} else {
EVCs <- EVC_local %>%
sf::st_filter(y = y, .predicate = st_intersects)
}
cols_to_g <- colnames(y %>% sf::st_drop_geometry())
# Summarise the resulting data
EVCs_summarised <- sf::st_intersection(y, EVCs) %>%
dplyr::mutate(AREASQM = sf::st_area(.)) %>%
dplyr::group_by_at(dplyr::vars(!!!cols_to_g, !!!group_vars)) %>%
dplyr::summarise(AREASQM = sum(AREASQM, na.rm = TRUE)) %>%
dplyr::ungroup() %>%
dplyr::mutate(evc_perc = as.numeric((100*AREASQM)/4e6))# Add empty EVC spot for missing EVC dat
options(vicmap.chunk_limit = o)
return(EVCs_summarised)
}
# fire severity
fire_severity <- terra::rast(here::here("_posts", "2021-09-21-southernarkfire", "data", "2020_bushfire", "ANZVI0803008638_FireSeverityFinal_20200414.tif"))
gbc <- gips_bound %>% sf::st_transform(3111) %>% vect()
fire_severity_c <- mask(crop(fire_severity, gbc), gbc)
values(fire_severity_c)[values(fire_severity_c) == 0] <- NA
Load Data
If Data is updated on Proofsafe, or the transects are updated you can change the write_data object to TRUE and provide it with new URLs. Otherwise data will be read from ARI’s postgresql database. Additionally, we also write data for EVCs and ALA presence-only records, this is for convenience and avoids having to pull down this spatial data from external databases everytime. We also have sections of code for processing earlier collected data from Greater Glider surveys (contact: Jemma Cripps, ARI).
Show code
if(write_data) {
raw_headers <- readr::read_csv("http://proofsafe.com.au/_xdata/1815-416-1159-Header-36fbe40a173a4327bc397dcd7bfc5202.csv")
raw_animals <- readr::read_csv("http://proofsafe.com.au/_xdata/1815-416-1161-Animals-e842915dc4e74d2f9af23f9c2aed6eee.csv")
DBI::dbWriteTable(con, Id(schema = "koala",
table = "raw_headers"),
raw_headers,
overwrite = TRUE)
DBI::dbWriteTable(con, Id(schema = "koala",
table = "raw_animals"),
raw_animals,
overwrite = TRUE)
# sf::st_write(koala_sites, "koala_sites.kml", delete_layer = TRUE)
koala_transects <- sf::read_sf(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "koala_transects.kml")) %>%
sf::st_transform(4283) %>%
tidyr::separate(col = Name, into = c("Site", "Transect"))
st_write(obj = koala_transects,
dsn = con,
layer = Id(schema = "koala", table = "koala_transects"),
delete_layer = TRUE)
}
if(write_evcs) {
evcs <- list()
for(i in 1:nrow(koala_transects)) {
evcs[[i]] <- extract_evc.vicmap_promise(y = koala_transects[i, ])
}
names(evcs) <- paste(koala_transects$Site, koala_transects$Transect, sep = ".")
evc_flat <- lapply(evcs, function(x) {
x %>%
group_by(Site, Transect) %>%
summarise(BIOREGION = paste(unique(BIOREGION), collapse = " | "),
EVC = paste(unique(X_GROUPNAME), collapse = " | "))
})
koala_evc_df <- bind_rows(evc_flat)
st_write(obj = koala_evc_df,
dsn = con,
layer = Id(schema = "koala", table = "koala_evc_df"),
delete_layer = TRUE)
}
if(ala_write) {
koala_ala_records <- galah::ala_occurrences(taxa = select_taxa("Phascolarctos cinereus"), mint_doi = TRUE,
filters = select_filters(stateProvince = "Victoria",
year = seq(2012, 2022))) %>%
filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>%
sf::st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4283)
st_write(obj = koala_ala_records,
dsn = con,
layer = Id(schema = "koala", table = "koala_ala_records"),
delete_layer = TRUE)
}
if(gg_nocturnal_write) {
ARI_GG_Sites <- sf::st_read(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data","GGNocturnal", "ARI_RFA_GG_spotlight_sites"))
ARI_GG_FIRE_sites <- sf::st_read(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data","GGNocturnal", "GG_PostFire_IPA_BBRR_RBR_transect_pts_all"))
ARI_Bogies <- sf::st_read(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data","GGNocturnal", "ARI_Greater_Glider_transects_Strathbogies_ALL"))
ARI_BOA <- sf::st_read(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data","GGNocturnal", "BoA_sites_surveyed"))
GG_proofsafe_header <- readr::read_csv("http://proofsafe.com.au/_xdata/1815-135-338-Header-5842d1220412472aa336a03f2b4bb785.csv")
GG_proofsafe_animals <- readr::read_csv("http://proofsafe.com.au/_xdata/1815-135-340-Animals-bfb42396299e4c54b02ac011932c2c96.csv")
ARI_BOA_form <- ARI_BOA %>%
mutate(transect_length = as.numeric(sf::st_length(geometry))) %>%
dplyr::select(SiteID = SiteID,
transect_length) %>%
sf::st_centroid()
ARI_Bogies_filt <- ARI_Bogies %>%
mutate(transect_length = as.numeric(sf::st_length(geometry)),
id = as.character(id)) %>%
dplyr::select(SiteID = id,
transect_length) %>%
sf::st_centroid()
ARI_GG_Sites_filt <- ARI_GG_Sites %>%
filter(StartEnd == "START") %>%
mutate(SiteID = case_when(Site == 0 ~ stringr::str_remove(PtName, "END|START"),
round(Site) == Site ~ as.character(as.integer(Site)),
round(Site) != Site ~ as.character(Site))) %>%
dplyr::select(SiteID,
transect_length = TransectLe)
ARI_GG_FIRE_sites_dist <- ARI_GG_FIRE_sites %>%
filter(start_end %in% c("start", "end")) %>%
group_by(project_si) %>%
summarise(transect_length = st_distance(geometry[1], geometry[2])) %>%
ungroup() %>%
sf::st_drop_geometry()
proofsafe_fire <- GG_proofsafe_header %>%
filter(stringr::str_detect(string = H0_SiteID, pattern = "BBRR|IPA|G1-|G2-|G3-|G4-")) %>%
dplyr::select(SiteID = H0_SiteID, H0_Transect, H1_Notes, Data_File_Id, H0_Observer1)
ARI_GG_FIRE_sites_formatted <- ARI_GG_FIRE_sites %>%
filter(start_end == "mid") %>%
left_join(ARI_GG_FIRE_sites_dist, by = "project_si") %>%
mutate(SiteID = case_when(project == "BBRR" ~ project_si,
project == "IPA" ~ project_si,
project == "RBR" ~ site),
SiteID = toupper(SiteID)) %>%
dplyr::select(SiteID, transect_length)
proofsafe_fire_formatted <- proofsafe_fire %>%
mutate(SiteID = stringr::str_replace(SiteID, "IPA_|IPA1_|IPA2_|IPA3_|IPA4_|IPA5_", "IPA-"),
SiteID = stringr::str_replace(SiteID, "_", "-"),
SiteID = toupper(SiteID),
SiteID = stringr::str_replace(SiteID, "BBRR-F", "BBRR-"),
SiteID = stringr::str_replace(SiteID, "IPA-299.1", "IPA-299"))
fire_joined <- left_join(proofsafe_fire_formatted, ARI_GG_FIRE_sites_formatted) %>%
mutate(transect_length = as.numeric(transect_length))
proofsafe_all <- GG_proofsafe_header %>%
dplyr::select(SiteID = H0_SiteID, H0_Transect, H1_Notes, Data_File_Id, H0_Observer1)
ARI_GG_Sites_filt2 <- ARI_GG_Sites_filt %>%
mutate(transect_length = as.numeric(transect_length)) %>%
bind_rows(ARI_BOA_form) %>%
bind_rows(ARI_Bogies_filt)
all_joined <- left_join(proofsafe_all, ARI_GG_Sites_filt) %>%
mutate(transect_length = as.numeric(transect_length)) %>%
filter(!is.na(transect_length)) %>%
bind_rows(fire_joined)
GG_proofsafe_animals_formatted <- GG_proofsafe_animals %>%
dplyr::mutate(Species = dplyr::case_when(I0_Species == "Other" ~ I0_Animal_sp_other,
TRUE ~ I0_Species),
TreeSpecies = dplyr::case_when(I0_Tree_species == "Other" ~ as.character(I0_Tree_sp_other),
TRUE ~ I0_Tree_species)) %>%
dplyr::select(Data_File_Id,
Animal_number = I0_Animal,
Species,
ObservationTime = I0_AnimObsTime,
ObservationType = I0_SeenHeard,
Direction = I0_L_or_R_of_trans,
Waypoint = I0_Waypoint_no_,
Distance = I0_Distance_to_animal,
Bearing = I0_Bearing_to_A_,
DistanceAlongTransect = I0_Dist_F_Transect,
TreeSpecies,
SeenByOther = I0_SeenX2,
VBAUpload = I0_VBAUse,
Comments = I0_Comments
) %>%
filter(Species == "Koala" & !is.na(Distance) & Distance <= 25) %>%
mutate(SeenByOther = case_when(is.na(SeenByOther) ~ "no",
TRUE ~ tolower(SeenByOther)))
GG_Gippsland_Survey_Data <- all_joined %>%
sf::st_as_sf() %>%
left_join(GG_proofsafe_animals_formatted, by = "Data_File_Id") %>%
group_by(SiteID, transect_length) %>%
summarise(Both = sum(SeenByOther == "yes")/2,
Obs1 = sum(H0_Observer1 == 1 & SeenByOther == "no"),
Obs2 = sum(H0_Observer1 == 1 & SeenByOther == "no")) %>%
ungroup() %>%
replace_na(list(Both = 0, Obs1 = 0, Obs2 = 0)) %>%
left_join(all_joined %>%
dplyr::select(SiteID) %>%
unique()) %>%
mutate(`Plot_area` = as.numeric(0.05*(transect_length/1000)*100),
Habitat = "Native") %>%
sf::st_as_sf() %>%
sf::st_cast("POINT") %>%
mutate(Method = "Nocturnal Double Count") %>%
sf::st_transform(3111)
GG_Gippsland_Survey_Data <- GG_Gippsland_Survey_Data[!st_is_empty(GG_Gippsland_Survey_Data),]
# GG_Gippsland_Survey_Data_Thinned <- spThin::thin(loc.data = GG_Gippsland_Survey_Data %>%
# st_coordinates() %>%
# as.data.frame() %>%
# mutate(SPEC = "Koala"),
# write.log.file = FALSE,
# lat.col = "Y",
# long.col = "X",
# write.files = FALSE,
# spec.col = "SPEC",
# thin.par = 1, # 1 km is same resolution as the environment covariates
# reps = 10,
# locs.thinned.list.return = TRUE)
st_write(obj = GG_Gippsland_Survey_Data,
dsn = con,
layer = Id(schema = "koala", table = "GG_Gippsland_Survey_Data"),
delete_layer = TRUE)
}
if(koala_translocations_write) {
koala_translocations <- sf::read_sf(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "historic_koala_translocations.geojson")) %>%
mutate(value = 1)
KoalaPredLayersAll <- terra::rast(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "KoalaPredLayers.tif"))
translocations_raster <- terra::rasterize(vect(koala_translocations %>%
sf::st_transform(3111)),
KoalaPredLayersAll[[1]],
fun = "max", na.rm = TRUE, background = NA, overwrite = TRUE)
translocation_dist_rast <- terra::distance(translocations_raster, grid = TRUE) %>%
mask(KoalaPredLayersAll[[1]])
names(translocation_dist_rast) <- "DistanceToTranslocation"
}
# Koala Sites
koala_sites <- sf::read_sf(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "Koala_Monitoring_Sites_editedDW")) %>%
sf::st_transform(4283) %>%
dplyr::select(name = SIteID)
raw_headers <- dplyr::tbl(con, dbplyr::in_schema("koala", "raw_headers")) %>%
collect()
koala_transects <- sf::st_read(con, layer = Id(schema = "koala", table = "koala_transects"))
all_transects <- koala_data_transect_view(raw_headers = raw_headers, koala_transects = koala_transects) %>%
mutate(EffortTime = EndTime - StartTime,
SiteTransect = paste(Site, Transect, sep = "."))
raw_animals <- dplyr::tbl(con, dbplyr::in_schema("koala", "raw_animals")) %>%
collect()
all_koala_observations <- koala_data_presence_view(raw_headers, raw_animals, koala_transects) %>%
filter(Species == "Koala") %>%
mutate(EffortTime = EndTime - StartTime,
OtherWaypoint = case_when(SeenByOther == "Yes" & User_Name %in% c("Justin", "Hayley") ~ gsub( " .*$", "", Comments)),
SiteTransect = paste(Site, Transect, sep = "."),
UniqueObs = case_when(Method == "Double Observer" & VBAUpload == "yes" ~ TRUE,
Method == "Double Observer" & VBAUpload != "yes" ~ FALSE,
Method == "Thermal" ~ TRUE))
# Data simplify cut
# koala_export <- list()
#
# koala_export[["Site_Information"]] <- all_transects %>%
# filter(Site < 5 & Method == "Double Observer", ObserverPosition == 1) %>%
# dplyr::select(-geometry, - Data_File_Id, -User_Name, -Observer, -ObserverPosition, -EffortTime, -SiteTransect, -Description)
#
# koala_export[["Koala_Observations"]] <- all_koala_observations %>%
# filter(Site < 5 & Method == "Double Observer") %>%
# dplyr::select(-Data_File_Id, -User_Name, -Observer, -Temperature, -EndTime, -StartTime, -EffortTime, -VBAUpload, -SiteTransect, -geometry, -UniqueObs, -Notes) %>%
# mutate(KoalaID = case_when(ObserverPosition == 1 ~ paste0("Koala_", Waypoint),
# ObserverPosition == 2 & SeenByOther == "Yes" ~ paste0("Koala_", OtherWaypoint),
# TRUE ~ paste0("Koala_", Waypoint))) %>%
# dplyr::select(KoalaID, everything())
#
#
# library(openxlsx)
# # create workbook
# wb <- createWorkbook()
#
# #Iterate the same way as PavoDive, slightly different (creating an anonymous function inside Map())
# Map(function(data, nameofsheet){
#
# addWorksheet(wb, nameofsheet)
# writeData(wb, nameofsheet, data)
#
# }, koala_export, names(koala_export))
#
# ## Save workbook to excel file
# saveWorkbook(wb, file = "data/KoalaExampleExport.xlsx", overwrite = TRUE)
Data Summary
Survey Effort
We surveyed a total of 70 transects across 39 sites. Three methods were used to survey: Double Observer Distance sampling, single observer thermal-assisted ground surveys and single-observer spotlight search. The survey effort (in both time and distance) varied between methods.
Show code
translength = all_transects %>%
dplyr::select(Method, SiteTransect, geometry) %>%
distinct() %>%
sf::st_as_sf() %>%
group_by(Method) %>%
summarise(`Total Transect Length` = sf::st_length(geometry) %>% sum() %>% units::set_units(km) %>% round(2)) %>%
sf::st_drop_geometry()
at <- koala_transects %>%
group_by(Site) %>%
summarise(l = sf::st_length(geometry) %>%
sum() %>%
units::set_units(m) %>% round(0))
atmin <- min(at$l)
atmax <- max(at$l)
atmean <- round(mean(at$l))
survey_summary <- all_transects %>%
group_by(Method) %>%
summarise(Sites = length(unique(Site)),
Transects = length(unique(SiteTransect)),
`Time Effort` = lubridate::seconds_to_period(sum(EffortTime)),
`Distance Effort` = sf::st_length(geometry) %>% sum() %>% set_units(km) %>% round(2),
`Average Speed` = (drop_units(`Distance Effort`)/(lubridate::time_length(`Time Effort`, "hours"))) %>% set_units("km/h") %>% round(2)) %>%
ungroup() %>%
left_join(translength, by = "Method")
survey_summary %>%
kbl(caption = "Summary of survey effort by method") %>%
kable_styling()
| Method | Sites | Transects | Time Effort | Distance Effort | Average Speed | Total Transect Length |
|---|---|---|---|---|---|---|
| Double Observer | 39 | 69 | 1d 16H 5M 53S | 63.89 [km] | 1.59 [km/h] | 31.94 [km] |
| Spotlight | 14 | 27 | 7H 42M 0S | 11.91 [km] | 1.55 [km/h] | 11.61 [km] |
| Thermal | 15 | 28 | 1d 3H 28M 41S | 11.75 [km] | 0.43 [km/h] | 11.75 [km] |
The mean combined transect distance at sites was 836, with the minimum at 443 and the maximum at 1679
Koala Detections
Koalas were detected at only 5 out of 39 sites. The detection of koalas using the double observer method and thermal method are not directly comparable as they were not conducted immediately after one another. Thus the assumption of static animals is violated. The double-observer method allows for the calculation of detection probability based on the variable detection of individual koalas because transects were walked ~15 minutes apart and assumes koalas to remain stationary in this time.
Show code
nokoalas <- anti_join(all_transects, all_koala_observations, by = c("Method", "SiteTransect")) %>%
group_by(Method, SiteTransect) %>%
summarise(`Koalas` = 0L,
`Koalas Seen By Both` = 0L)
koala_per_sites <- all_koala_observations %>%
group_by(Method, SiteTransect) %>%
summarise(`Koalas` = sum(UniqueObs == TRUE),
`Koalas Seen By Both` = as.integer(sum(SeenByOther == "Yes")/2)) %>%
ungroup() %>%
mutate(`Koalas Seen By Both` = case_when(Method != "Double Observer" ~ NA_integer_,
TRUE ~ `Koalas Seen By Both`)) %>%
bind_rows(nokoalas) %>%
as.data.table() %>%
dcast(SiteTransect ~ Method, value.var = c("Koalas", "Koalas Seen By Both")) %>%
dplyr::select(`Site - Transect` = SiteTransect,
`Individual Koalas Seen By Double Observer Method` = `Koalas_Double Observer`,
`Individual Koalas Seen By Both Observers` = `Koalas Seen By Both_Double Observer`,
`Individual Koalas Seen By Thermal Method` = `Koalas_Thermal`) %>%
arrange(desc(`Individual Koalas Seen By Double Observer Method`), desc(`Individual Koalas Seen By Thermal Method`)) %>%
add_row(cbind(tibble(`Site - Transect` = "Total"), t(as.data.frame(colSums(.[,2:4], na.rm = T))) %>% `rownames<-`(NULL)))
koala_per_sites %>%
kbl(caption = "Koalas observed from Double Observer and Thermal Surveys") %>%
kable_styling() %>%
add_header_above(c(" " = 1, "Double Observer Method" = 2, "Thermal Method" = 1)) %>%
row_spec(nrow(koala_per_sites), bold = TRUE)
| Site - Transect | Individual Koalas Seen By Double Observer Method | Individual Koalas Seen By Both Observers | Individual Koalas Seen By Thermal Method |
|---|---|---|---|
| 3.2 | 7 | 4 | 12 |
| 3.1 | 6 | 1 | 13 |
| 4.2 | 4 | 2 | 6 |
| 4.1 | 2 | 2 | 5 |
| 55.1 | 1 | 0 | 2 |
| 1.1 | 1 | 0 | 0 |
| 1.3 | 1 | 0 | 0 |
| 1.2 | 0 | 0 | 0 |
| 11.1 | 0 | 0 | 0 |
| 11.2 | 0 | 0 | 0 |
| 12.1 | 0 | 0 | 0 |
| 14.1 | 0 | 0 | 0 |
| 14.2 | 0 | 0 | 0 |
| 2.1 | 0 | 0 | 0 |
| 2.2 | 0 | 0 | 0 |
| 20.1 | 0 | 0 | 0 |
| 23.1 | 0 | 0 | 0 |
| 23.2 | 0 | 0 | 0 |
| 28.2 | 0 | 0 | 0 |
| 28.3 | 0 | 0 | 0 |
| 42.1 | 0 | 0 | 0 |
| 44.1 | 0 | 0 | 0 |
| 45.1 | 0 | 0 | 0 |
| 45.2 | 0 | 0 | 0 |
| 53.1 | 0 | 0 | 0 |
| 53.3 | 0 | 0 | 0 |
| 55.2 | 0 | 0 | 0 |
| 10.1 | 0 | 0 | NA |
| 10.2 | 0 | 0 | NA |
| 13.1 | 0 | 0 | NA |
| 13.2 | 0 | 0 | NA |
| 18.1 | 0 | 0 | NA |
| 21.1 | 0 | 0 | NA |
| 21.2 | 0 | 0 | NA |
| 21.3 | 0 | 0 | NA |
| 22.1 | 0 | 0 | NA |
| 22.2 | 0 | 0 | NA |
| 22.3 | 0 | 0 | NA |
| 25.1 | 0 | 0 | NA |
| 26.1 | 0 | 0 | NA |
| 26.2 | 0 | 0 | NA |
| 28.1 | 0 | 0 | NA |
| 32.1 | 0 | 0 | NA |
| 32.2 | 0 | 0 | NA |
| 34.1 | 0 | 0 | NA |
| 40.1 | 0 | 0 | NA |
| 40.2 | 0 | 0 | NA |
| 41.1 | 0 | 0 | NA |
| 41.2 | 0 | 0 | NA |
| 43.1 | 0 | 0 | NA |
| 49.1 | 0 | 0 | NA |
| 5.1 | 0 | 0 | NA |
| 5.2 | 0 | 0 | NA |
| 50.1 | 0 | 0 | NA |
| 51.1 | 0 | 0 | NA |
| 53.4 | 0 | 0 | NA |
| 54.1 | 0 | 0 | NA |
| 54.2 | 0 | 0 | NA |
| 56.1 | 0 | 0 | NA |
| 56.2 | 0 | 0 | NA |
| 57.1 | 0 | 0 | NA |
| 58.1 | 0 | 0 | NA |
| 58.2 | 0 | 0 | NA |
| 6.1 | 0 | 0 | NA |
| 7.1 | 0 | 0 | NA |
| 7.2 | 0 | 0 | NA |
| 8.1 | 0 | 0 | NA |
| 8.2 | 0 | 0 | NA |
| 9.1 | 0 | 0 | NA |
| 53.2 | NA | NA | 1 |
| Total | 22 | 9 | 39 |
Koalas were detected at varying distances from transects, and this varied between survey methods. With thermal-assisted surveys showing detections at further distances. Unfortunately, for this study we do not have a large enough sample size to accurately use double-observer distance sampling to generate detection functions.
Show code
all_koala_observations %>%
ggplot() +
geom_histogram(aes(x = HorizontalDistance), bins = 10, fill = "Grey50") +
facet_wrap(~Method, scales = "fixed", nrow = 2) +
labs(x = "Horizontal Distance to Koala", y = "Number of Koalas Observed") +
theme_bw()
Despite the lack of records, there is some suggestive evidence that detectibility may vary depending on distance from observer and height of animal. Overall, there is variability in detection based on horizontal and vertical distance.
Show code
all_koala_observations %>%
filter(Method == "Double Observer" & VBAUpload == "yes") %>%
ggplot() +
geom_point(aes(x = HorizontalDistance, y = Height, fill = SeenByOther, size = SeenByOther), shape = 21) +
labs(x = "Horizontal Distance", fill = "Seen By Both Observers", size = "Seen By Both Observers") +
theme_bw()
Figure 1: Double observer distance sampling varies over horizontal and vertical distances
Map of surveys and detections
Koalas were detected at sites around Mallacoota and Gelantipy. With two sites at Gelantipy having very high koala densities. Surveys largely took place within the ‘East Gippsland Lowlands’, ‘East Gippsland Uplands’ and ‘Monaro Tablelands’ bioregions.
Show code
koala_records <-readRDS(paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/historic_koala_ala_records.rds"))) %>%
sf::st_intersection(gips_bound) %>%
sf::st_coordinates() %>%
as.data.frame() %>%
mutate(Presence = "Historical Presence")
koala_map <- koala_per_sites %>% left_join(koala_transects %>%
mutate(`Site - Transect` = paste(Site, Transect, sep = ".")),
by = c("Site - Transect")) %>%
group_by(Site) %>%
summarise(Koala = coalesce(sum(c(`Individual Koalas Seen By Double Observer Method`)#,
# `Individual Koalas Seen By Thermal Method`)
, na.rm = TRUE), 0),
geometry = sf::st_centroid(sf::st_union(geometry))) %>%
ungroup() %>%
mutate(Presence = case_when(Koala != 0 ~ "Present",
TRUE ~ "Absent")) %>%
sf::st_as_sf() %>%
cbind(sf::st_coordinates(.)) %>%
rename(latitude = Y, longitude = X)
koala_map$Presence <- factor(koala_map$Presence, levels = c("Absent", "Present"))
koala_map <- arrange(koala_map, Presence)
koala_map %>%
ggplot() +
geom_sf(data = gips_bound) +
geom_sf(data = evc_sub, aes(fill = BIOREGION), alpha = 0.5) +
geom_point(aes(x = X, y = Y, colour = Presence), shape = 22, size = 2, alpha = 1, data = koala_records, stroke = 1, fill = "white") +
geom_point(aes(x = longitude, y = latitude, colour = Presence), shape = 19, size = 2) +
scale_color_manual(values = c(Absent = "Black", Present = "Grey95", `Historical Presence` = "Black")) +
guides(color = guide_legend(override.aes = list(shape = c(19,19,22), fill = "white") ) ) +
theme_void()
Figure 2: Current surveys and historical presence compared against historical records (since 2007) and the various bioregions in East Gippsland
Surveys were based on an existing model of koala density. However this model did not encorperate any surveys in East Gippsland and thus relied upon the habitat suitability (e.g. preferred eucalypt feed trees). Most surveys were restricted to this existing layer, however due to the likely inaccuracy of the extent of non-zero koala density from the original model, some sites outside this layer were sampled if habitat or historical records suggested koalas may be present (e.g. Mallacoota).
Show code
koala_pred <- terra::rast(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "Abundance_prediction_native.tif"))
terra::values(koala_pred)[terra::values(koala_pred) == 0] <- NA
koala_pred_crop <- terra::mask(koala_pred, vect(gips_bound %>%
sf::st_transform(3111))) %>%
crop(gips_bound %>% sf::st_transform(3111))
# mapview(koala_pred, na.color = "transparent")
koala_map %>%
ggplot() +
geom_sf(data = gips_bound, fill = "White") +
geom_spatraster(data = koala_pred_crop, aes(fill = Abundance_prediction_native)) +
scale_fill_gradient2(na.value = "transparent", low = delwp_cols[["Environment"]], mid = delwp_cols[["Teal"]], high = delwp_cols[["Navy"]], midpoint = 0.5, limits = c(0,1), guide = guide_colourbar(title = "Predicted Koala \nDensity (ha)\n")) +
# scale_fill_gradientn(
# colors = hcl.colors(10, "RdBu", rev = TRUE),
# na.value = NA, guide = guide_colourbar(title = "Predicted Koala \nDensity\n")) +
geom_point(aes(x = X, y = Y, colour = Presence), shape = 22, size = 2, alpha = 1, data = koala_records, stroke = 1, fill = "white") +
geom_point(aes(x = longitude, y = latitude, colour = Presence, shape = Presence), size = 3) +
scale_color_manual(values = c(Absent = "Black", Present = "Red", `Historical Presence` = "Black")) +
guides(color = guide_legend(override.aes = list(shape = c(19,17,22), fill = "white")), shape = FALSE) +
theme_void()
Figure 3: Current surveys and historical presence compared against historical records (since 2007) and the existing predicted modelled koala density
Show code
fireimpact <- extract(fire_severity_c, vect(koala_transects %>% st_transform(3111)), method = "simple") %>%
group_by(ID) %>%
summarise(`Transect Fire Impact` = paste(unique(ANZVI0803008638_FireSeverityFinal_20200414), collapse = ", "),
`Maximum Fire Impact` = max(ANZVI0803008638_FireSeverityFinal_20200414, na.rm = TRUE))
fireimpact$Site <- as.integer(koala_transects$Site)
fireimpact$Transect <- koala_transects$Transect
fireimpact_sum <- fireimpact %>%
group_by(Site) %>%
summarise(`Maximum Impact` = max(`Maximum Fire Impact`, na.rm = T)) %>%
ungroup() %>%
mutate(`Maximum Impact` = case_when(is.infinite(`Maximum Impact`) ~ 0,
TRUE ~ `Maximum Impact`))
nofire <- sum(fireimpact_sum$`Maximum Impact` == 0)
lowfire <- sum(fireimpact_sum$`Maximum Impact` == 3)
midfire <- sum(fireimpact_sum$`Maximum Impact` == 4)
highfire <- sum(fireimpact_sum$`Maximum Impact` >= 5)
fireimpact_sum %>%
arrange(Site) %>%
dplyr::select(Site, `Maximum Impact`) %>%
kbl(caption = "Maximum Fire Severity at Each Site") %>%
kable_styling("striped", full_width = FALSE) %>%
collapse_rows(1) %>%
scroll_box(height = "1000px", width = "800px")
| Site | Maximum Impact |
|---|---|
| 1 | 3 |
| 2 | 5 |
| 3 | 5 |
| 4 | 0 |
| 5 | 4 |
| 6 | 4 |
| 7 | 0 |
| 8 | 0 |
| 9 | 5 |
| 10 | 4 |
| 11 | 4 |
| 12 | 4 |
| 13 | 0 |
| 14 | 3 |
| 18 | 6 |
| 20 | 0 |
| 21 | 0 |
| 22 | 3 |
| 23 | 4 |
| 25 | 0 |
| 26 | 5 |
| 28 | 0 |
| 32 | 0 |
| 34 | 0 |
| 40 | 0 |
| 41 | 0 |
| 42 | 5 |
| 43 | 3 |
| 44 | 0 |
| 45 | 5 |
| 49 | 4 |
| 50 | 6 |
| 51 | 6 |
| 53 | 5 |
| 54 | 5 |
| 55 | 5 |
| 56 | 0 |
| 57 | 0 |
| 58 | 0 |
In total, 16 sites had no fire impact, 4 sites had low fire impact, 7 sites had medium fire impact and 12 sites had high fire impact (Department of Environment 2022).
Vegetation
Based on the Eucalyptus habitat assessment we compare the predicted presence of seven eucalyptus koala food trees E. camaldulensis, E. camphora, E. cypellocarpa, E. globulus, E. ovata, E. rubida and E. viminalis. The predicted percentage is the proportion of the 1km-square grid that is covered by that species (e.g. hectares). False positive error is defined as having a predicted proportion above 5% (5ha) and not being detected by surveyers.
Show code
euc_stub <- here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "Eucs")
e.camald <- rast(paste0(euc_stub, '/E_camaldulensis_1000_m.tif'))
e.camph <- rast(paste0(euc_stub, '/E_camphora_1000_m.tif'))
e.cyp <- rast(paste0(euc_stub, '/E_cypellocarpa_1000_m.tif'))
e.glob <- rast(paste0(euc_stub, '/E_globulus_1000_m.tif'))
e.ova <- rast(paste0(euc_stub, '/E_ovata_1000_m.tif'))
e.rub <- rast(paste0(euc_stub, '/E_rubida_1000_m.tif'))
e.vim <- rast(paste0(euc_stub, '/E_viminalis_1000_m.tif'))
euc_stack <- c(e.camald,
e.camph,
e.cyp,
e.glob,
e.ova,
e.rub,
e.vim)
names(euc_stack) <- c("E. camaldulensis",
"E. camphora",
"E. cypellocarpa",
"E. globulus",
"E. ovata",
"E. rubida",
"E. viminalis")
eucalypt_detections <- readxl::read_excel(paste0(euc_stub, "/EastGippslandTreeID.xlsx")) %>%
filter(!(Site %in% c(34,45)))
eucalypt_detections_long <- as.data.table(eucalypt_detections %>%
dplyr::select(-Order,
-Date,
-Description)) %>%
melt(id.vars = c("Site", "Dominant")) %>%
mutate(Site = as.character(Site))
# for each site get predicted raw koala layers
transects_site <- koala_transects %>%
filter(!(Site %in% c("34","45"))) %>%
group_by(Site) %>%
summarise() %>%
sf::st_transform(3111)
euc_prop <- bind_cols(transects_site[,"Site"] %>% st_drop_geometry(),
extract(euc_stack, vect(transects_site)) %>%
group_by(ID) %>%
summarise_all("sum")) %>%
dplyr::select(-ID) %>%
as.data.table() %>%
melt(id.vars = "Site", value.name = "Predicted Hectares")
euc_prop_combined <- left_join(euc_prop,
eucalypt_detections_long,
by = c("Site", "variable")) %>%
rename(Species = variable)
euc_prop_combined_summary <- euc_prop_combined %>%
group_by(Site, Species, Dominant) %>%
summarise(`Predicted Percentage` = mean(`Predicted Hectares`),
`Ground Truthed` = sum(as.numeric(value), na.rm = TRUE)) %>%
ungroup() %>%
mutate(Site = as.integer(Site)) %>%
arrange(Site) %>%
mutate(Error = case_when(`Predicted Percentage` > 5 & `Ground Truthed` == 0 ~ "False Positive",
`Predicted Percentage` == 0 & `Ground Truthed` == 1 ~ "False Negative",
TRUE ~ "No Error"))
euc_prop_combined_summary %>%
kbl(caption = "Modelled Eucalypt Feed Trees Ground Truthing. For thw 7 eucalypt food trees we compare the predicted presence against a ground truth assessment of presence") %>%
kable_styling(bootstrap_options = "striped") %>%
collapse_rows(1:2) %>%
scroll_box(height = "1000px", width = "800px")
| Site | Species | Dominant | Predicted Percentage | Ground Truthed | Error |
|---|---|---|---|---|---|
| 1 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. globulus | 1.84 | 0 | No Error | |
| E. cypellocarpa | E. globulus | 69.87 | 1 | No Error | |
| E. globulus | E. globulus | 17.58 | 1 | No Error | |
| E. ovata | E. globulus | 3.16 | 0 | No Error | |
| E. rubida | E. globulus | 8.66 | 0 | False Positive | |
| E. viminalis | E. globulus | 27.10 | 0 | False Positive | |
| 2 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. globoidea | 1.10 | 0 | No Error | |
| E. cypellocarpa | E. globoidea | 75.08 | 1 | No Error | |
| E. globulus | E. globoidea | 14.16 | 1 | No Error | |
| E. ovata | E. globoidea | 2.02 | 0 | No Error | |
| E. rubida | E. globoidea | 5.90 | 0 | False Positive | |
| E. viminalis | E. globoidea | 22.00 | 0 | False Positive | |
| 3 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. viminalis | 2.34 | 1 | No Error | |
| E. cypellocarpa | E. viminalis | 30.99 | 0 | False Positive | |
| E. globulus | E. viminalis | 8.41 | 0 | False Positive | |
| E. ovata | E. viminalis | 2.69 | 0 | No Error | |
| E. rubida | E. viminalis | 7.85 | 0 | False Positive | |
| E. viminalis | E. viminalis | 21.82 | 1 | No Error | |
| 4 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. viminalis | 0.04 | 0 | No Error | |
| E. cypellocarpa | E. viminalis | 0.18 | 0 | No Error | |
| E. globulus | E. viminalis | 0.08 | 0 | No Error | |
| E. ovata | E. viminalis | 0.02 | 1 | No Error | |
| E. rubida | E. viminalis | 0.43 | 0 | No Error | |
| E. viminalis | E. viminalis | 0.35 | 1 | No Error | |
| 5 | E. camaldulensis | NA | 0.09 | 0 | No Error |
| E. camphora | MIX | 0.44 | 0 | No Error | |
| E. cypellocarpa | MIX | 16.09 | 1 | No Error | |
| E. globulus | MIX | 1.98 | 0 | No Error | |
| E. ovata | MIX | 5.05 | 0 | False Positive | |
| E. rubida | MIX | 0.11 | 0 | No Error | |
| E. viminalis | MIX | 5.35 | 0 | False Positive | |
| 6 | E. camaldulensis | NA | 0.02 | 0 | No Error |
| E. camphora | MIX | 0.12 | 0 | No Error | |
| E. cypellocarpa | MIX | 20.87 | 0 | False Positive | |
| E. globulus | MIX | 1.89 | 0 | No Error | |
| E. ovata | MIX | 3.22 | 0 | No Error | |
| E. rubida | MIX | 0.15 | 0 | No Error | |
| E. viminalis | MIX | 2.37 | 0 | No Error | |
| 7 | E. camaldulensis | NA | 1.34 | 0 | No Error |
| E. camphora | E. globoidea | 0.43 | 0 | No Error | |
| E. cypellocarpa | E. globoidea | 44.26 | 0 | False Positive | |
| E. globulus | E. globoidea | 3.41 | 0 | No Error | |
| E. ovata | E. globoidea | 13.69 | 0 | False Positive | |
| E. rubida | E. globoidea | 1.00 | 0 | No Error | |
| E. viminalis | E. globoidea | 19.89 | 0 | False Positive | |
| 8 | E. camaldulensis | NA | 0.91 | 0 | No Error |
| E. camphora | MIX | 0.04 | 0 | No Error | |
| E. cypellocarpa | MIX | 25.70 | 0 | False Positive | |
| E. globulus | MIX | 2.79 | 0 | No Error | |
| E. ovata | MIX | 9.15 | 0 | False Positive | |
| E. rubida | MIX | 0.18 | 0 | No Error | |
| E. viminalis | MIX | 14.12 | 0 | False Positive | |
| 9 | E. camaldulensis | NA | 0.11 | 0 | No Error |
| E. camphora | E. botryoides | 0.00 | 0 | No Error | |
| E. cypellocarpa | E. botryoides | 22.12 | 0 | False Positive | |
| E. globulus | E. botryoides | 2.27 | 0 | No Error | |
| E. ovata | E. botryoides | 5.43 | 0 | False Positive | |
| E. rubida | E. botryoides | 0.02 | 0 | No Error | |
| E. viminalis | E. botryoides | 5.64 | 0 | False Positive | |
| 10 | E. camaldulensis | NA | 0.13 | 0 | No Error |
| E. camphora | E. botryoides | 0.00 | 0 | No Error | |
| E. cypellocarpa | E. botryoides | 19.91 | 0 | False Positive | |
| E. globulus | E. botryoides | 1.19 | 0 | No Error | |
| E. ovata | E. botryoides | 16.34 | 0 | False Positive | |
| E. rubida | E. botryoides | 0.00 | 0 | No Error | |
| E. viminalis | E. botryoides | 7.20 | 0 | False Positive | |
| 11 | E. camaldulensis | NA | 0.05 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 20.60 | 0 | False Positive | |
| E. globulus | MIX | 1.24 | 0 | No Error | |
| E. ovata | MIX | 4.24 | 0 | No Error | |
| E. rubida | MIX | 0.02 | 0 | No Error | |
| E. viminalis | MIX | 1.45 | 0 | No Error | |
| 12 | E. camaldulensis | NA | 0.12 | 0 | No Error |
| E. camphora | MIX | 0.01 | 0 | No Error | |
| E. cypellocarpa | MIX | 16.60 | 0 | False Positive | |
| E. globulus | MIX | 1.92 | 0 | No Error | |
| E. ovata | MIX | 4.05 | 0 | No Error | |
| E. rubida | MIX | 0.03 | 0 | No Error | |
| E. viminalis | MIX | 3.05 | 0 | No Error | |
| 13 | E. camaldulensis | NA | 3.07 | 0 | No Error |
| E. camphora | E. globoidea | 0.32 | 0 | No Error | |
| E. cypellocarpa | E. globoidea | 53.96 | 0 | False Positive | |
| E. globulus | E. globoidea | 10.81 | 0 | False Positive | |
| E. ovata | E. globoidea | 10.07 | 0 | False Positive | |
| E. rubida | E. globoidea | 0.48 | 0 | No Error | |
| E. viminalis | E. globoidea | 20.72 | 0 | False Positive | |
| 14 | E. camaldulensis | NA | 1.55 | 0 | No Error |
| E. camphora | MIX | 0.35 | 0 | No Error | |
| E. cypellocarpa | MIX | 75.33 | 0 | False Positive | |
| E. globulus | MIX | 7.72 | 1 | No Error | |
| E. ovata | MIX | 11.22 | 0 | False Positive | |
| E. rubida | MIX | 0.86 | 0 | No Error | |
| E. viminalis | MIX | 19.89 | 0 | False Positive | |
| 18 | E. camaldulensis | NA | 0.58 | 0 | No Error |
| E. camphora | E. polyanthemos | 1.12 | 0 | No Error | |
| E. cypellocarpa | E. polyanthemos | 59.55 | 0 | False Positive | |
| E. globulus | E. polyanthemos | 18.89 | 0 | False Positive | |
| E. ovata | E. polyanthemos | 3.30 | 0 | No Error | |
| E. rubida | E. polyanthemos | 4.36 | 0 | No Error | |
| E. viminalis | E. polyanthemos | 17.48 | 0 | False Positive | |
| 20 | E. camaldulensis | NA | 1.73 | 0 | No Error |
| E. camphora | E. globoidea | 0.45 | 0 | No Error | |
| E. cypellocarpa | E. globoidea | 22.37 | 1 | No Error | |
| E. globulus | E. globoidea | 1.48 | 0 | No Error | |
| E. ovata | E. globoidea | 8.98 | 0 | False Positive | |
| E. rubida | E. globoidea | 0.71 | 0 | No Error | |
| E. viminalis | E. globoidea | 12.49 | 0 | False Positive | |
| 21 | E. camaldulensis | NA | 2.38 | 0 | No Error |
| E. camphora | E. polyanthemos | 0.15 | 0 | No Error | |
| E. cypellocarpa | E. polyanthemos | 12.45 | 1 | No Error | |
| E. globulus | E. polyanthemos | 3.12 | 0 | No Error | |
| E. ovata | E. polyanthemos | 5.63 | 0 | False Positive | |
| E. rubida | E. polyanthemos | 1.27 | 0 | No Error | |
| E. viminalis | E. polyanthemos | 12.92 | 0 | False Positive | |
| 22 | E. camaldulensis | NA | 1.44 | 0 | No Error |
| E. camphora | E. cypellocarpa | 0.08 | 0 | No Error | |
| E. cypellocarpa | E. cypellocarpa | 14.10 | 1 | No Error | |
| E. globulus | E. cypellocarpa | 2.91 | 0 | No Error | |
| E. ovata | E. cypellocarpa | 4.29 | 0 | No Error | |
| E. rubida | E. cypellocarpa | 0.69 | 0 | No Error | |
| E. viminalis | E. cypellocarpa | 9.82 | 0 | False Positive | |
| 23 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. sieberi | 0.07 | 0 | No Error | |
| E. cypellocarpa | E. sieberi | 67.47 | 0 | False Positive | |
| E. globulus | E. sieberi | 11.84 | 0 | False Positive | |
| E. ovata | E. sieberi | 1.26 | 0 | No Error | |
| E. rubida | E. sieberi | 0.23 | 0 | No Error | |
| E. viminalis | E. sieberi | 6.38 | 0 | False Positive | |
| 25 | E. camaldulensis | NA | 4.73 | 0 | No Error |
| E. camphora | E. polyanthemos | 1.19 | 0 | No Error | |
| E. cypellocarpa | E. polyanthemos | 29.39 | 1 | No Error | |
| E. globulus | E. polyanthemos | 7.62 | 0 | False Positive | |
| E. ovata | E. polyanthemos | 9.65 | 0 | False Positive | |
| E. rubida | E. polyanthemos | 2.26 | 0 | No Error | |
| E. viminalis | E. polyanthemos | 22.87 | 0 | False Positive | |
| 26 | E. camaldulensis | NA | 0.77 | 0 | No Error |
| E. camphora | E. polyanthemos | 0.08 | 0 | No Error | |
| E. cypellocarpa | E. polyanthemos | 22.75 | 1 | No Error | |
| E. globulus | E. polyanthemos | 3.69 | 0 | No Error | |
| E. ovata | E. polyanthemos | 2.70 | 0 | No Error | |
| E. rubida | E. polyanthemos | 0.46 | 0 | No Error | |
| E. viminalis | E. polyanthemos | 6.29 | 0 | False Positive | |
| 28 | E. camaldulensis | NA | 1.25 | 0 | No Error |
| E. camphora | E. globoidea | 0.32 | 0 | No Error | |
| E. cypellocarpa | E. globoidea | 34.99 | 1 | No Error | |
| E. globulus | E. globoidea | 3.85 | 0 | No Error | |
| E. ovata | E. globoidea | 12.87 | 0 | False Positive | |
| E. rubida | E. globoidea | 0.43 | 0 | No Error | |
| E. viminalis | E. globoidea | 18.29 | 0 | False Positive | |
| 32 | E. camaldulensis | NA | 0.04 | 0 | No Error |
| E. camphora | MIX | 15.61 | 0 | False Positive | |
| E. cypellocarpa | MIX | 46.23 | 0 | False Positive | |
| E. globulus | MIX | 7.37 | 0 | False Positive | |
| E. ovata | MIX | 2.14 | 0 | No Error | |
| E. rubida | MIX | 56.44 | 1 | No Error | |
| E. viminalis | MIX | 79.23 | 0 | False Positive | |
| 40 | E. camaldulensis | NA | 1.81 | 0 | No Error |
| E. camphora | E. botryoides | 0.14 | 0 | No Error | |
| E. cypellocarpa | E. botryoides | 24.96 | 0 | False Positive | |
| E. globulus | E. botryoides | 4.37 | 0 | No Error | |
| E. ovata | E. botryoides | 10.55 | 0 | False Positive | |
| E. rubida | E. botryoides | 0.15 | 0 | No Error | |
| E. viminalis | E. botryoides | 17.17 | 0 | False Positive | |
| 41 | E. camaldulensis | NA | 2.48 | 0 | No Error |
| E. camphora | MIX | 0.33 | 0 | No Error | |
| E. cypellocarpa | MIX | 54.66 | 1 | No Error | |
| E. globulus | MIX | 4.35 | 0 | No Error | |
| E. ovata | MIX | 17.90 | 0 | False Positive | |
| E. rubida | MIX | 0.67 | 0 | No Error | |
| E. viminalis | MIX | 26.11 | 0 | False Positive | |
| 42 | E. camaldulensis | NA | 0.38 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 19.92 | 0 | False Positive | |
| E. globulus | MIX | 1.87 | 0 | No Error | |
| E. ovata | MIX | 9.28 | 0 | False Positive | |
| E. rubida | MIX | 0.08 | 0 | No Error | |
| E. viminalis | MIX | 7.18 | 0 | False Positive | |
| 43 | E. camaldulensis | NA | 0.51 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 8.08 | 0 | False Positive | |
| E. globulus | MIX | 0.93 | 0 | No Error | |
| E. ovata | MIX | 9.58 | 0 | False Positive | |
| E. rubida | MIX | 0.02 | 0 | No Error | |
| E. viminalis | MIX | 16.02 | 0 | False Positive | |
| 44 | E. camaldulensis | NA | 1.13 | 0 | No Error |
| E. camphora | MIX | 0.05 | 0 | No Error | |
| E. cypellocarpa | MIX | 42.45 | 1 | No Error | |
| E. globulus | MIX | 4.14 | 0 | No Error | |
| E. ovata | MIX | 10.90 | 0 | False Positive | |
| E. rubida | MIX | 0.18 | 0 | No Error | |
| E. viminalis | MIX | 14.26 | 0 | False Positive | |
| 49 | E. camaldulensis | NA | 1.19 | 0 | No Error |
| E. camphora | E. polyanthemos | 0.21 | 0 | No Error | |
| E. cypellocarpa | E. polyanthemos | 35.00 | 1 | No Error | |
| E. globulus | E. polyanthemos | 6.70 | 0 | False Positive | |
| E. ovata | E. polyanthemos | 3.43 | 0 | No Error | |
| E. rubida | E. polyanthemos | 0.38 | 0 | No Error | |
| E. viminalis | E. polyanthemos | 10.15 | 0 | False Positive | |
| 50 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 15.71 | 0 | False Positive | |
| E. globulus | MIX | 1.13 | 0 | No Error | |
| E. ovata | MIX | 4.89 | 0 | No Error | |
| E. rubida | MIX | 0.00 | 0 | No Error | |
| E. viminalis | MIX | 1.28 | 0 | No Error | |
| 51 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 30.40 | 0 | False Positive | |
| E. globulus | MIX | 2.48 | 0 | No Error | |
| E. ovata | MIX | 6.99 | 0 | False Positive | |
| E. rubida | MIX | 0.01 | 0 | No Error | |
| E. viminalis | MIX | 1.59 | 0 | No Error | |
| 53 | E. camaldulensis | NA | 0.05 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 15.63 | 0 | False Positive | |
| E. globulus | MIX | 1.39 | 0 | No Error | |
| E. ovata | MIX | 10.64 | 0 | False Positive | |
| E. rubida | MIX | 0.00 | 0 | No Error | |
| E. viminalis | MIX | 3.21 | 0 | No Error | |
| 54 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 23.76 | 0 | False Positive | |
| E. globulus | MIX | 1.96 | 0 | No Error | |
| E. ovata | MIX | 4.80 | 0 | No Error | |
| E. rubida | MIX | 0.00 | 0 | No Error | |
| E. viminalis | MIX | 2.27 | 0 | No Error | |
| 55 | E. camaldulensis | NA | 0.21 | 0 | No Error |
| E. camphora | MIX | 0.00 | 0 | No Error | |
| E. cypellocarpa | MIX | 15.52 | 1 | No Error | |
| E. globulus | MIX | 1.59 | 0 | No Error | |
| E. ovata | MIX | 5.57 | 0 | False Positive | |
| E. rubida | MIX | 0.00 | 0 | No Error | |
| E. viminalis | MIX | 4.86 | 0 | No Error | |
| 56 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. croajingolensis | 8.79 | 0 | False Positive | |
| E. cypellocarpa | E. croajingolensis | 31.43 | 0 | False Positive | |
| E. globulus | E. croajingolensis | 4.95 | 0 | No Error | |
| E. ovata | E. croajingolensis | 0.62 | 0 | No Error | |
| E. rubida | E. croajingolensis | 38.40 | 1 | No Error | |
| E. viminalis | E. croajingolensis | 43.97 | 0 | False Positive | |
| 57 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. pauciflora | 5.24 | 0 | False Positive | |
| E. cypellocarpa | E. pauciflora | 17.20 | 0 | False Positive | |
| E. globulus | E. pauciflora | 5.35 | 0 | False Positive | |
| E. ovata | E. pauciflora | 0.70 | 0 | No Error | |
| E. rubida | E. pauciflora | 18.81 | 0 | False Positive | |
| E. viminalis | E. pauciflora | 21.05 | 1 | No Error | |
| 58 | E. camaldulensis | NA | 0.00 | 0 | No Error |
| E. camphora | E. mannifera | 6.68 | 0 | False Positive | |
| E. cypellocarpa | E. mannifera | 30.84 | 0 | False Positive | |
| E. globulus | E. mannifera | 16.77 | 0 | False Positive | |
| E. ovata | E. mannifera | 1.04 | 0 | No Error | |
| E. rubida | E. mannifera | 37.63 | 0 | False Positive | |
| E. viminalis | E. mannifera | 35.99 | 0 | False Positive |
Show code
euc_prop_combined_summary2 <- euc_prop_combined_summary %>%
group_by(Error) %>%
summarise(Sites = length(unique(Site)),
Count = n())
euc_summary3 <- euc_prop_combined_summary %>%
group_by(Error, Species) %>%
summarise(Sites = length(unique(Site)),
`Average Prediction` = mean(`Predicted Percentage`)) %>%
as.data.table() %>%
dcast(Species ~ Error, value.var = c("Sites", "Average Prediction")) %>%
mutate(`Sites_False Positive` = tidyr::replace_na(`Sites_False Positive`, 0),
`Average Prediction_False Positive` = tidyr::replace_na(`Average Prediction_False Positive` , 0)) %>%
arrange(desc(`Sites_False Positive`)) %>%
mutate(`Error Rate` = `Sites_False Positive`/(`Sites_No Error`+ `Sites_False Positive`)) %>%
dplyr::select(Species,
`False Positive` = `Sites_False Positive`,
`No Error` = `Sites_No Error`,
`Error Rate`,
`Average False Positive Predicted Hectares` = `Average Prediction_False Positive`)
euc_summary3 %>%
kbl(caption = "Modelled Eucalypt Feed Trees Ground Truthing. For the 7 Eucalyptus food trees we assess the number sites where the type of error occured", digits = 2) %>%
kable_styling(bootstrap_options = "striped") %>%
collapse_rows(1:2)
| Species | False Positive | No Error | Error Rate | Average False Positive Predicted Hectares |
|---|---|---|---|---|
| E. viminalis | 26 | 11 | 0.70 | 19.17 |
| E. cypellocarpa | 23 | 14 | 0.62 | 31.37 |
| E. ovata | 19 | 18 | 0.51 | 9.97 |
| E. globulus | 9 | 28 | 0.24 | 10.42 |
| E. rubida | 5 | 32 | 0.14 | 15.77 |
| E. camphora | 4 | 33 | 0.11 | 9.08 |
| E. camaldulensis | 0 | 37 | 0.00 | 0.00 |
Based on the assessment of modelled food trees and those detected with ground surveys we can see that all species apart from E. camaldulensis have error rates above 10 %. E. viminalis and E. cypellocarpa appear to be especially problematic with many sites showing predictions that could not be ground truthed. The high amount of false positives and the lack of false negatives tell us one or both of two things:
- Modelled Eucalyptus koala food tree species is overestimated.
- Detection rates for the koala food tree species during ground-truthing is low.
Koala Model
The koala density model consists of a multi-step process to (i) determine habitat suitability and (ii) predict density within the suitable habitat area.
Habitat suitability threshold
A previous koala density model was written by Geoff Heard and Dave Ramsey. Here we update the model with new data from this survey as well as other presence-only records. In comparison to the previously written Maxent model here we utilise a more extensive model selection process as well as including more recent detection data. Using the SDMtune package (Vignali et al. 2020) we undertake the following model selection and refinement process:
- Download presence-only data from atlas of living australia (using the
galahpackage) (Stevenson, Westgate, and Newman 2022). Data is restricted to Victoria and since 2007.
- Spatially-thin data (1km buffer) using
spThin(Aiello-Lammens et al. 2015) to prevent oversampling in incidental detection hotspots, and have the same resolution as environmental predictors.
- Select background points randomly across Victoria. The number of background points is set to two times the number of presence detections.
- Prepare a ‘Sample With Data’ object using the presence, background and environmental raster.
- Generate a ‘base model’ that includes all environmental variables and a regularisation parameter of 2.
- Cut data into training (0.8) and testing portions (0.2).
- Select variables (using
varSel()) for inclusion in the model based on removing correlated variables (pearson correlation > 0.8).
- Reduce variables for low influence variables (below 10%) that do not improve AUC.
- Tune hyperparameters of regularisation (reg) and feature classes (fc) by using a
randomSearch(), select top model based on AICc.
- Select top performing models from tuned model selection.
- From the predicted top model generate a threshold model that equalizes
test sensitivity and specificity.
With regard to the Eucalyptus layers the preference layers include:
High preference: Eucalyptus viminalis
Medium preference: Eucalyptus ovata + Eucalyptus globulus + Eucalyptus camaldulensis
Low preference: Eucalyptus rubida + Eucalyptus camphora + Eucalyptus cypellocarpa
Show code
predictionfile <- paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/KoalaMaxentmap_terra.tif")
KoalaPredLayers <- raster::stack(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "KoalaPredLayers.tif"))
KoalaPredLayers <- subset(KoalaPredLayers, 1:21)
# KoalaPredLayers[[22]] <- raster(sqrt(translocation_dist_rast))
area <- vic_bound %>%
sf::st_union() %>%
sf::st_transform(3111) %>%
sf::st_buffer(25000) %>%
sf::st_simplify(dTolerance = 10000) %>%
sf::st_transform(4283) %>%
sf::st_as_sf()
if(!file.exists(predictionfile)) {
# Read in precompiled environmental data. In order to generate this stack see:
# source()
historic_koala_ala_records <- galah_call() %>%
galah_identify("Phascolarctos cinereus") %>%
galah::galah_geolocate(area) %>%
galah_filter(year >= 2007) %>%
atlas_occurrences(mint_doi = TRUE) %>%
filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>%
sf::st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4283) %>%
mutate(eventDate = as.Date(eventDate))
historic_koala_ala_records$eventDate <- as.Date(historic_koala_ala_records$eventDate)
historic_koala_ala_records <- historic_koala_ala_records %>%
filter(eventDate < as.Date("2022-06-30"))
koala_ala_points <- sf::st_coordinates(historic_koala_ala_records %>% sf::st_transform(3111))
koala_thinned_points <- spThin::thin(loc.data = koala_ala_points %>%
as.data.frame() %>%
mutate(SPEC = "Koala"),
lat.col = "Y",
long.col = "X",
write.files = FALSE,
spec.col = "SPEC",
write.log.file = FALSE,
thin.par = 1, # 1 km is same resolution as the environment covariates
reps = 10,
locs.thinned.list.return = TRUE)
koala_thinned_points_top <- koala_thinned_points[[length(koala_thinned_points)]]
# Select background points
koala_background <- randomPoints(KoalaPredLayers,
n = nrow(koala_thinned_points_top)*2,
p = koala_thinned_points_top,
ext = raster::extent(KoalaPredLayers),
tryf = 10)
# Generate data for SDMtune model
koala_maxent_data <- prepareSWD(species = "Koala",
env = KoalaPredLayers,
p = koala_thinned_points_top,
a = koala_background)
# Conduct 5-fold CV
folds <- randomFolds(koala_maxent_data, k = 5,
only_presence = TRUE,
seed = 25)
# Default parameters except regularisation increased to 3
base_model <- train(method = "Maxent",
data = koala_maxent_data,
folds = folds, reg = 2)
base_modelNF <- train(method = "Maxent",
data = koala_maxent_data, reg = 2)
cat("Training AUC: ", auc(base_model))
cat("Testing AUC: ", auc(base_model, test = TRUE))
#### Select variables ####
train_test <- trainValTest(koala_maxent_data,
test = 0.2,
only_presence = TRUE,
seed = 25)
train <- train_test[[1]]
test <- train_test[[2]]
set.seed(25)
bg_p <- randomPoints(KoalaPredLayers,
n = 500,
ext = raster::extent(KoalaPredLayers))
bg <- prepareSWD(species = "Bgs", a = bg_p, env = KoalaPredLayers)
koala_selected_variables_model <- varSel(base_modelNF,
metric = "auc",
test = test,
bg4cor = bg,
method = "pearson",
cor_th = 0.8,
permut = 5)
reduced_variables_model <- reduceVar(koala_selected_variables_model,
th = 10,
metric = "auc",
test = test,
permut = 5,
use_jk = TRUE)
h <- list(reg = seq(2, 5, 0.5),
fc = c("l", "lq", "lh", "lp", "lqp", "lqph"))
koala_rstune <- randomSearch(reduced_variables_model,
hypers = h,
metric = "aic",
test = test,
env = KoalaPredLayers,
pop = 10,
seed = 25)
TopKoalaModel <- train("Maxent",
data = koala_rstune@models[[1]]@data,
fc = koala_rstune@results[1, 1],
reg = koala_rstune@results[1, 2],
folds = folds)
TopKoalaModelNoCV <- train("Maxent",
data = koala_rstune@models[[1]]@data,
fc = koala_rstune@results[1, 1],
reg = koala_rstune@results[1, 2])
TopKoalaModelVi <- maxentVarImp(TopKoalaModel)
TopKoalaModelVi %>%
kbl("html", caption = "Variable importance of covariates included in the top model explaining Koala Distribution", digits = 3) %>%
kable_styling(full_width = FALSE)
cont_plots <- TopKoalaModel@data@data %>% names()
plots <- lapply(cont_plots, function(x) {
plotResponse(TopKoalaModel,
var = x,
type = "cloglog",
only_presence = TRUE,
marginal = FALSE,
color = "black") +
ylab("Habitat suitability")
})
cowplot::plot_grid(plotlist = plots)
stack_filtered <- KoalaPredLayers[[TopKoalaModel@data@data %>% names()]]
KoalaMaxentmap <- predict(TopKoalaModelNoCV,
data = stack_filtered,
type = "cloglog",
progress = "text")
KoalaMaxentmap_terra <- terra::rast(KoalaMaxentmap) %>%
terra::mask(terra::rast(fasterize::fasterize(vic_bound, raster = KoalaMaxentmap)))
writeRaster(KoalaMaxentmap_terra,
file = predictionfile,
overwrite = TRUE)
saveRDS(test, paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/maxenttest.rds")))
saveRDS(koala_thinned_points_top, paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/koala_thinned_points_top.rds")))
saveRDS(TopKoalaModelNoCV, paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/TopKoalaModelNoCV.rds")))
saveRDS(historic_koala_ala_records, paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/historic_koala_ala_records.rds")))
} else {
KoalaMaxentmap_terra <- terra::rast(predictionfile)
koala_thinned_points_top <- readRDS(paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/koala_thinned_points_top.rds")))
test <- readRDS(paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/maxenttest.rds")))
TopKoalaModelNoCV <- readRDS(paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/TopKoalaModelNoCV.rds")))
historic_koala_ala_records <- readRDS(paste0(paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/historic_koala_ala_records.rds")))
}
statewideMaxentPred <- ggplot() +
geom_spatraster(data = KoalaMaxentmap_terra, aes(fill = layer), na.rm = TRUE) +
geom_sf(data = vic_bound, fill = "grey", inherit.aes = FALSE, alpha = 0, colour = "grey70") +
# geom_sf(data = regions_3111, fill = "grey", inherit.aes = FALSE, alpha = 0.2, colour = "grey70") +
# geom_sf(data = statewide_water_simp, colour = "#2c7fb8", inherit.aes = FALSE, alpha = 0.5) +
# geom_tile(aes(fill = value), na.rm = TRUE) +
scale_fill_gradient2(na.value = "transparent", low = delwp_cols[["Environment"]], mid = delwp_cols[["Teal"]], high = delwp_cols[["Navy"]], midpoint = 0.5, limits = c(0,1)) +
ylab("Latitude") +
xlab("Longitude") +
labs(fill = "Habitat \nSuitability") +
theme(axis.title = element_text(size = 40),
axis.text = element_text(size = 30),
legend.title = element_text(size=40),
legend.text = element_text(size=30),
legend.key.size = unit(6, 'cm')) +
theme_bw()
ths <- thresholds(TopKoalaModelNoCV, type = "cloglog", test = test)
sampledTH <- terra::spatSample(KoalaMaxentmap_terra, prod(dim(KoalaMaxentmap_terra)), "regular", as.raster = TRUE)
sampledTH_dat <- as.data.frame(sampledTH, xy=TRUE, na.rm = FALSE)
sampledTH_dat$mask <- factor(cut(sampledTH_dat$layer,
breaks=c(0,ths[5, 2], 1), labels = c(NA_integer_, 1)))
sampledTH_dat$TH <- as.factor(cut(sampledTH_dat$layer,
breaks=c(0,ths[5, 2], 1), labels = c(0, 1)))
thresholdRaster <- terra::rast(raster::rasterFromXYZ(sampledTH_dat[,c(1:5)] %>% rename(HabitatSuitability = layer), res = res(sampledTH), crs = crs(sampledTH)))
values(thresholdRaster$TH) <- factor(values(thresholdRaster$TH))
# Turn 15 years of records to raster
koalarecords_mask <- fasterize::fasterize(historic_koala_ala_records %>%
sf::st_transform(3111) %>%
mutate(Pres = 1) %>%
st_buffer(25000),
raster = raster(thresholdRaster$TH), field = "Pres", background = NA) %>%
rast() %>%
terra::mask(x = ., mask = thresholdRaster$mask)
add(thresholdRaster) <- koalarecords_mask
names(thresholdRaster) <- c(names(thresholdRaster)[1:3], "koalarecords_mask")
values(thresholdRaster$koalarecords_mask) <- factor(values(thresholdRaster$koalarecords_mask))
dfvals <- data.frame(to = 1:length(values(thresholdRaster$mask)), mask = values(thresholdRaster$mask))
ad <- adjacent(thresholdRaster$mask,
which(!is.na(values(thresholdRaster$mask))),
directions=16,
pairs=TRUE) %>%
as.data.frame() %>%
left_join(dfvals, by = "to") %>%
filter(!is.na(mask)) %>%
group_by(from) %>%
summarise(cells = n())
#### KDE Estimate of distribution ####
st_kde <- function(points,cellsize, bandwith, extent = NULL){
if(is.null(extent)){
extent_vec <- st_bbox(points)[c(1,3,2,4)]
} else{
extent_vec <- st_bbox(extent)[c(1,3,2,4)]
}
n_y <- ceiling((extent_vec[4]-extent_vec[3])/cellsize)
n_x <- ceiling((extent_vec[2]-extent_vec[1])/cellsize)
extent_vec[2] <- extent_vec[1]+(n_x*cellsize)-cellsize
extent_vec[4] <- extent_vec[3]+(n_y*cellsize)-cellsize
coords <- st_coordinates(points)
matrix <- MASS::kde2d(coords[,1],coords[,2],h = bandwith,n = c(n_x,n_y),lims = extent_vec)
raster::raster(matrix)
}
hist_3111 <- historic_koala_ala_records %>%
sf::st_transform(3111)
# bw.x <- stats::bw.SJ(koala_thinned_points_top[,1],
# nb = 1000,
# method = "dpi")
bw.x.nrd <- bandwidth.nrd(hist_3111 %>% st_coordinates() %>% .[,1])
# bw.y <- stats::bw.SJ(koala_thinned_points_top[,2],
# nb = 1000,
# method = "dpi")
bw.y.nrd <- bandwidth.nrd(hist_3111 %>% st_coordinates() %>% .[,2])
points_kde <- st_kde(hist_3111,
cellsize = 1000,
bandwith = bandwidth.nrd(hist_3111 %>% st_coordinates()), extent = terra::ext(thresholdRaster$mask))
pred_vals <- extract(points_kde, hist_3111)
q995 <- quantile(pred_vals, 0.005, na.rm = T)
points_kde_cut <- points_kde
values(points_kde_cut)[values(points_kde_cut) >= q995] <- as.factor("1")
values(points_kde_cut)[values(points_kde_cut) < q995] <- NA
points_kde_cut <- rast(points_kde_cut)
terra::crs(points_kde_cut) <- terra::crs(thresholdRaster$mask)
points_kde_cut <- resample(points_kde_cut, thresholdRaster$mask)
points_kde_cut <- mask(points_kde_cut, thresholdRaster$mask)
names(points_kde_cut) <- "KDE Distribution"
add(thresholdRaster) <- points_kde_cut
values(thresholdRaster$`KDE Distribution`) <- factor(values(thresholdRaster$`KDE Distribution`), levels = c(NA, "1"))
thresholdfile <- paste0(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data"), "/KoalaMaxentTH.tif")
writeRaster(thresholdRaster,
file = thresholdfile,
overwrite = TRUE)
#
# th_poly <- terra::as.polygons(thresholdRaster, trunc = TRUE, dissolve = TRUE) %>%
# sf::st_as_sf()
#
# cc_poly <- concaveman::concaveman(points = th_poly %>% filter(cut == 1) %>% sf::st_cast("MULTIPOINT"), concavity = 0.9)
#
# threshold_smooth <- smoothr::smooth(th_poly, method = "ksmooth", smoothness = 2)raster
historic_recs_crop <- historic_koala_ala_records %>%
sf::st_transform(3111) %>%
sf::st_intersection(vic_bound) %>%
sf::st_coordinates() %>%
as.data.frame()
statewideMaxentPredTH <- ggplot() +
# geom_sf(data = regions_3111, fill = "grey", inherit.aes = FALSE, alpha = 0.2, colour = "grey70") +
# geom_sf(data = statewide_water_simp, colour = "#2c7fb8", inherit.aes = FALSE, alpha = 0.5) +
geom_spatraster(aes(fill = TH),
na.rm = TRUE, data = thresholdRaster[["TH"]]) +
geom_point(data = historic_recs_crop,
aes(x = X, y = Y), shape = 21, size = 1, fill = "red") +
scale_fill_manual(na.value = "transparent", values = c(`1` = unname(delwp_cols["Navy"]), `0` = unname(delwp_cols["Environment"]))) +
ylab("Latitude") +
xlab("Longitude") +
labs(fill = "Habitat \nSuitability") +
theme(axis.title = element_text(size = 40),
axis.text = element_text(size = 30),
legend.title = element_text(size=40),
legend.text = element_text(size=30),
legend.key.size = unit(6, 'cm')) +
theme_bw()
statewideMaxentPredTH_recs <- ggplot() +
geom_sf(data = vic_bound, fill = unname(delwp_cols["Environment"]), inherit.aes = FALSE, alpha = 1, colour = "grey70") +
# geom_sf(data = statewide_water_simp, colour = "#2c7fb8", inherit.aes = FALSE, alpha = 0.5) +
geom_spatraster(aes(fill = `KDE Distribution`),
na.rm = TRUE, data = thresholdRaster[["KDE Distribution"]]) +
# geom_point(data = historic_recs_crop,
# aes(x = X, y = Y), shape = 21, size = 1, fill = "red") +
scale_fill_manual(na.value = "transparent", values = c(`1` = unname(delwp_cols["Navy"]), `0` = unname(delwp_cols["Environment"]))) +
ylab("Latitude") +
xlab("Longitude") +
labs(fill = "Habitat \nSuitability") +
theme(axis.title = element_text(size = 40),
axis.text = element_text(size = 30),
legend.title = element_text(size=40),
legend.text = element_text(size=30),
legend.key.size = unit(6, 'cm')) +
theme_bw()
Our tuned top model has an AUC of 0.925. This improves on previous maxent models, which had a maximum AUC of 0.787. The total area across the state is 4,199,900 ha.
Based on the Maxent Habitat Suitability Model we can generate statewide predictions for habitat suitability. It should be noted that despite spatial-thinning issues remain regarding likely unequal sampling effort. Models of habitat suitability tend to perform better for Western Victoria distributions. Populations in Eastern Victoria; where fewer data points, and arguably fewer surveys have been conducted prove challanging for modelling. Several possible options on improving habitat models exist:
- Improve background point selection to account for bias in sampling.
- Include habitat and data for southern New South Wales, as these populations are likely to be similar to the Eastern/North-Eastern Victorian population.
- Addition of a wider variety of preference trees (e.g. yellow gum and iron bark).
Show code
statewideMaxentPred
Figure 4: Statewide habitat suitability for Koalas from a tuned Maxent model. Predictions are formed from koala presence records since 2007
Show code
cowplot::plot_grid(statewideMaxentPredTH, statewideMaxentPredTH_recs, nrow = 2, ncol = 1, labels = "AUTO")
Figure 5: Statewide habitat suitability threshold for Koalas from a tuned Maxent model. The threshold is based on equal test sensitivity and specificity. (A) shows the unrestricted threshold for all of Victoria while (B) shows the threshold restricted to areas within the Kernel Density Estimate
The threshold for the top model is based on the maximum test sensitivity plus specificity. The Cloglog value for this threshold is 0.179, above which, areas are considered as ‘suitable habitat’. The fraction of the state included in this threshold is 0.163. The omission rate of test data (i.e, koala records not used in the model is 1.2 %). The total area of the threshold raster (restricted to kernel density area) is 2,815,877 ha.
Point-process density models
Data Formatting
To rerun the model we are limited by the previously extracted datasets and thus we have regenerated statewide rasters of the environmental variables of interest. Below we describe the data preparation steps used to format model data:
Generate statewide rasters for all variables used in the model, or have the possibility of being used in the model. This includes regions and transformed variables for the eucalyptus layers as well as whether the area is plantation or native forest. The resolution of this data is set at 1km2 and contains the following layers: Rainfall, RainSeasonality, MTWQ, MTWM, MTDR, AnnualTempRange, TemperatureSeasonality, Altitude, Slope, Roughness, AllEuc, AllEucNBR, HighPrefEuc, HighPrefEucNBR, MedPrefEuc, MedPrefEucNBR, LowPrefEuc, LowPrefEucNBR, PlantationEuc, PlantationEucNBR, Nitrogen, hab.vals, Region, NativeProp and PlantationProp.
Format all records in a standardised way and ensure they are spatially explicit (sf objects). The format of these data objects should have overlapping core data columns of:
- Location/geometry
- Method
- Counts (obs1, obs2, both etc)
- Survey area
We have done this for previous data used in the last density model and for new data we format it similarly.
- Extract covariate data for all survey locations from raster.
- Based on the methods, split the data and format it according to the existing data that feeds into the STAN model.
- Additional data such as number of native and plantation cells (and their environmental covariates) can be easily calculated from masking the statewide raster with the threshold dataset.
Show code
KoalaPredLayersAll <- terra::rast(here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "KoalaPredLayers.tif"))
KoalaPredLayersAll_scaled <- terra::scale(subset(KoalaPredLayersAll, 1:20))
add(KoalaPredLayersAll_scaled) <- KoalaPredLayersAll$hab.vals
add(KoalaPredLayersAll_scaled) <- KoalaPredLayersAll$Region
add(KoalaPredLayersAll_scaled) <- KoalaPredLayersAll$NativeProp
add(KoalaPredLayersAll_scaled) <- KoalaPredLayersAll$PlantationProp
abundance.formula<- as.formula(~ hab.vals + Rainfall + Altitude + MTWQ + Slope +
HighPrefEucNBR + MedPrefEucNBR +
LowPrefEucNBR + PlantationEucNBR +
Rainfall * HighPrefEucNBR +
Rainfall * MedPrefEucNBR +
Rainfall * LowPrefEucNBR +
Rainfall * PlantationEucNBR +
Altitude * HighPrefEucNBR +
Altitude * MedPrefEucNBR +
Altitude * LowPrefEucNBR +
Altitude * PlantationEucNBR +
MTWQ * HighPrefEucNBR +
MTWQ * MedPrefEucNBR +
MTWQ * LowPrefEucNBR +
MTWQ * PlantationEucNBR)
# Extract covariates for the survey data
#### Single Count Data ####
single_count_data <- terra::extract(KoalaPredLayersAll_scaled,
vect(compiled_single_counts), method = "simple", factors = FALSE, list = FALSE) %>% as.data.frame()
single_count_data$Habitat <- compiled_single_counts$Habitat
single_count_data <- single_count_data %>%
mutate(hab.vals = case_when(Habitat == "Native" ~ 2,
TRUE ~ 1))
single_count_matrix <- model.matrix(abundance.formula, data=single_count_data[,-1])
#### Double Count Diurnal ####
## Add in NEW DATA ##
nokoalas_dist <- anti_join(all_transects, all_koala_observations %>%
filter(Method == "Double Observer"), by = c("Method", "Site")) %>%
filter(Method == "Double Observer") %>%
group_by(Site) %>%
summarise(`Both` = 0L,
`Obs1` = 0L,
`Obs2` = 0L)
koala_per_sites <- all_koala_observations %>%
filter(Method == "Double Observer" & HorizontalDistance <= 25) %>%
group_by(Site) %>%
summarise(Both = as.integer(sum(SeenByOther == "Yes")/2),
Obs1 = as.integer(sum(ObserverPosition == 1)),
Obs2 = as.integer(sum(ObserverPosition == 2))) %>%
ungroup() %>%
bind_rows(nokoalas_dist) %>%
arrange(Site) %>%
left_join(koala_map %>%
mutate(Site = as.numeric(Site)) %>%
dplyr::select(Site), by = "Site") %>%
sf::st_as_sf() %>%
sf::st_transform(3111) %>%
mutate(Site_ID = as.character(Site),
Habitat = "Native") %>%
dplyr::select(-Site)
# transectLength
do.gipp.trans.l <- all_transects %>%
mutate(Site = as.character(Site)) %>%
dplyr::select(Site_ID = Site, geometry) %>%
distinct() %>%
sf::st_as_sf() %>%
group_by(Site_ID) %>%
summarise(`Plot_area` = as.numeric(0.05*(sf::st_length(geometry) %>% sum() %>% units::set_units(km)))*100) %>%
sf::st_drop_geometry()
# Combine with existing data #
double_count_diurnal <- compiled_double_counts %>%
filter(Method == "Diurnal Double Count") %>%
bind_rows(koala_per_sites %>% left_join(do.gipp.trans.l, by = "Site_ID"))
double_count_diurnal_data <- terra::extract(KoalaPredLayersAll_scaled,
vect(double_count_diurnal),
method = "simple",
factors = FALSE,
list = FALSE) %>%
as.data.frame()
double_count_diurnal_data$Habitat <- double_count_diurnal$Habitat
double_count_diurnal_data <- double_count_diurnal_data %>%
mutate(hab.vals = case_when(Habitat == "Native" ~ 2,
TRUE ~ 1))
diurnal_double_count_matrix <- model.matrix(abundance.formula, data=double_count_diurnal_data[,-1])
#### Double count nocturnal ####
GG_Gippsland_Survey_Data <- sf::st_read(con, layer = Id(schema = "koala", table = "GG_Gippsland_Survey_Data")) %>%
sf::st_transform(3111) %>%
dplyr::select(Site_ID = SiteID,
Plot_area,
Both,
Obs1,
Obs2,
Method,
Habitat)
double_count_nocturnal <- compiled_double_counts %>%
filter(Method == "Nocturnal Double Count") %>%
bind_rows(GG_Gippsland_Survey_Data)
double_count_nocturnal_data <- terra::extract(KoalaPredLayersAll_scaled,
vect(double_count_nocturnal),
method = "simple",
factors = FALSE,
list = FALSE) %>%
as.data.frame()
double_count_nocturnal_data$Habitat <- double_count_nocturnal$Habitat
double_count_nocturnal_data <- double_count_nocturnal_data %>%
mutate(hab.vals = case_when(Habitat == "Native" ~ 2,
TRUE ~ 1))
nocturnal_double_count_matrix <- model.matrix(abundance.formula, data=double_count_nocturnal_data[,-1])
#### Native background data ####
vic.regions.df <- vicmap_query("datavic:VMADMIN_DELWP_REGION") %>%
collect() %>%
sf::st_transform(3111) %>%
mutate(DELWP_REGION_CODE = as.integer(DELWP_REGION_CODE)) %>%
st_drop_geometry() %>%
dplyr::select(Region = DELWP_REGION,
DELWP_REGION_CODE) %>%
mutate(DELWP_REGION_CODE2 = DELWP_REGION_CODE-1)
KoalaPredLayersAll_scaled.clip <- terra::mask(KoalaPredLayersAll_scaled, mask = thresholdRaster$`KDE Distribution`)
background_df <- terra::as.data.frame(KoalaPredLayersAll_scaled.clip, na.rm = FALSE, cells = TRUE) %>%
left_join(vic.regions.df)
# Get native veg in df format
native_df <- background_df %>%
filter(!is.na(NativeProp) & NativeProp >= 5 & !is.na(DELWP_REGION_CODE2))
native_mat <- model.matrix(abundance.formula, data=native_df)
plantation_df <- background_df %>%
filter(!is.na(PlantationProp) & PlantationProp > 1 & !is.na(DELWP_REGION_CODE2))
plantation_mat <- model.matrix(abundance.formula, data=plantation_df)
# all_sites <-
Stan Model
The stan model accepts a list with formatted data objects (processed above). The data includes the single count, double count (diurnal) and double count (nocturnal). Data is also provided for model predictions to native forest and plantation areas, as well as summarised region-specific estimates.
Note: For this current iteration of the model we are not including the records from the thermal surveys. They are single counts, with a different method to the other data, thus they are not directly congruent with the existing data format. As we obtain further data for further thermal surveys we will likely include this data.
The model is a point-process model with a zero-inflated poisson distribution.
Show code
data <- list(J=3,
ncb=ncol(single_count_matrix),
nsc=nrow(compiled_single_counts),
ndcd=nrow(double_count_diurnal),
ndcn=nrow(double_count_nocturnal),
ysc=compiled_single_counts$Count,
ydcd=double_count_diurnal[,c("Both", "Obs1", "Obs2")] %>%
sf::st_drop_geometry() %>%
as.matrix(),
ydcn=double_count_nocturnal[,c("Both", "Obs1", "Obs2")] %>%
sf::st_drop_geometry() %>%
as.matrix(),
hab = diurnal_double_count_matrix[,"hab.vals"],
Xsc=single_count_matrix,
Xdcd=diurnal_double_count_matrix,
Xdcn=nocturnal_double_count_matrix,
obs_sc=compiled_single_counts$Number_of_observers,
area_sc=compiled_single_counts$Plot_area,
area_dcd=double_count_diurnal$Plot_area,
area_dcn=double_count_nocturnal$Plot_area,
ncells=nrow(native_mat),
pcells=nrow(plantation_mat),
area_n=native_df$NativeProp,
area_p=plantation_df$PlantationProp,
Xpn=native_mat,
Xpp=plantation_mat,
n_reg_id=native_df$DELWP_REGION_CODE2,
p_reg_id=plantation_df$DELWP_REGION_CODE2)
The below chunk presents the STAN model. The model is largely unchanged except for limiting/truncating the total koala estimates to 5 koalas per hectare. Estimates above this number are likely sourced from model predictions in unobserved parameter space and are thus unrealistic.
Show code
data {
int<lower=0> J; // number of multinomial cells
int<lower=0> ncb; // number of covariates on beta
int<lower=0> nsc; // number of single count cells
int<lower=0> ndcd; // number of double count cells diurnal
int<lower=0> ndcn; // number of double count cells noctural
array[nsc] int<lower=0> ysc; // single count data
array[ndcd, J] int<lower=0> ydcd; // double count data day
array[ndcn, J] int<lower=0> ydcn; // double count data night
array[ndcd] int<lower=1> hab; // habitat indicator 1: Plantation, 2: native
matrix[nsc, ncb] Xsc; // Abundance covariates single counts
matrix[ndcd, ncb] Xdcd; // abundance covariates - double counts day
matrix[ndcn, ncb] Xdcn; // abundance covariates - double counts night
array[nsc] int<lower=0> obs_sc; // No. observers single counts
vector[nsc] area_sc; // offset single counts
vector[ndcd] area_dcd; // offset double counts day
vector[ndcn] area_dcn; // offset double counts night
int<lower=0> ncells; // No. native cells
int<lower=0> pcells; // No. plantation cells
vector[ncells] area_n; // offset native cells
vector[pcells] area_p; // offset plantation cells
matrix[ncells, ncb] Xpn; // native cell covariates
matrix[pcells, ncb] Xpp; // plantation cell covariates
array[ncells] int n_reg_id; // region indicators native
array[pcells] int p_reg_id; // region indicators planttaion
}
parameters {
vector[ncb] beta;
matrix<lower=0, upper=1>[2,2] Pdc;
vector<lower=0, upper=1>[2] Pdn;
real<lower=0> sd_sc;
real<lower=0> sd_dcd;
real<lower=0> sd_dcn;
vector[nsc] sc_raw;
vector[ndcd] dcd_raw;
vector[ndcn] dcn_raw;
}
transformed parameters {
vector<lower=0>[nsc] lambda_sc;
vector<lower=0>[ndcd] lambda_dcd;
vector<lower=0>[ndcn] lambda_dcn;
vector[nsc] eps_sc;
vector[ndcd] eps_dcd;
vector[ndcn] eps_dcn;
vector<lower=0, upper=1>[nsc] psc;
matrix[ndcd, J] mud;
matrix[ndcn, J] mun;
real<lower=0, upper=1> Pave;
eps_sc = sd_sc * sc_raw;
eps_dcd = sd_dcd * dcd_raw;
eps_dcn = sd_dcn * dcn_raw;
lambda_sc = exp(Xsc * beta + eps_sc) .* area_sc;
lambda_dcd = exp(Xdcd * beta + eps_dcd) .* area_dcd;
lambda_dcn = exp(Xdcn * beta + eps_dcn) .* area_dcn;
Pave = mean(Pdc[1:2]);
for(j in 1:nsc) psc[j] = 1 - (1 - Pave)^obs_sc[j];
// multinomial cell probs for double counts
for(i in 1:ndcd) {
mud[i,1] = lambda_dcd[i] * Pdc[hab[i],1] * Pdc[hab[i],2];
mud[i,2] = lambda_dcd[i] * Pdc[hab[i],1] * (1-Pdc[hab[i],2]);
mud[i,3] = lambda_dcd[i] * (1-Pdc[hab[i],1]) * Pdc[hab[i],2];
}
// multinomial cell probs for double counts
for(i in 1:ndcn) {
mun[i,1] = lambda_dcn[i] * Pdn[1] * Pdn[2];
mun[i,2] = lambda_dcn[i] * Pdn[1] * (1-Pdn[2]);
mun[i,3] = lambda_dcn[i] * (1-Pdn[1]) * Pdn[2];
}
}
model {
beta ~ normal(0, 2);
to_vector(Pdc) ~ beta(2, 2);
Pdn ~ beta(2, 2);
sd_sc ~ student_t(4, 0, 1);
sd_dcd ~ student_t(4, 0, 1);
sd_dcn ~ student_t(4, 0, 1);
sc_raw ~ std_normal();
dcd_raw ~ std_normal();
dcn_raw ~ std_normal();
for(i in 1:ndcd)
target += poisson_lpmf(ydcd[i] | mud[i]);
for(i in 1:ndcn)
target += poisson_lpmf(ydcn[i] | mun[i]);
target += poisson_lpmf(ysc | lambda_sc .* psc);
}
generated quantities{
array[nsc] int<lower=0> ysc_rep;
array[ndcd, J] int<lower=0> ydcd_mat; // double count data day
array[ndcn, J] int<lower=0> ydcn_mat;
array[ndcd] int<lower=0> ydcd_rep;
array[ndcn] int<lower=0> ydcn_rep;
vector<lower=0>[ncells] Nn;
vector<lower=0>[pcells] Np;
vector<lower=0>[ncells] counter_n;
vector<lower=0>[pcells] counter_p;
real<lower=0> Nhat_n;
real<lower=0> Nhat_p;
vector<lower=0>[6] Nn_region;
vector<lower=0>[6] Np_region;
real<lower=0> Nhat;
ysc_rep = poisson_rng(lambda_sc .* psc);
for(j in 1:ndcd) {
ydcd_mat[j,1:J] = poisson_rng(mud[j,1:J]);
ydcd_rep[j] = sum(ydcd_mat[j]);
}
for(j in 1:ndcn) {
ydcn_mat[j,1:J] = poisson_rng(mun[j,1:J]);
ydcn_rep[j] = sum(ydcn_mat[j]);
}
for(i in 1:6) Nn_region[i] = 0;
for(i in 1:ncells) {
Nn[i] = poisson_log_rng(Xpn[i] * beta);
counter_n[i] = 0;
while(Nn[i] > 5) {
Nn[i] = poisson_log_rng(Xpn[i] * beta);
counter_n[i] += 1;
if(counter_n[i] > 100) break;
}
Nn_region[n_reg_id[i]] += Nn[i] * area_n[i];
}
for(i in 1:6) Np_region[i] = 0;
for(i in 1:pcells) {
Np[i] = poisson_log_rng(Xpp[i] * beta);
counter_p[i] = 0;
while(Np[i] > 5){
Np[i] = poisson_log_rng(Xpp[i] * beta);
counter_p[i] += 1;
if(counter_p[i] > 100) break;
}
Np_region[p_reg_id[i]] += Np[i] * area_p[i];
}
Nhat_n = sum(Nn .* area_n);
Nhat_p = sum(Np .* area_p);
Nhat = Nhat_n + Nhat_p;
}Show code
if(!file.exists("Outputs/Koala_model_outputs.stan")) {
# MCMC settings
ni <- 1000
nt <- 1
nb <- 1000
nc <- 4
inits <- lapply(1:nc, function(i) list(beta=rep(0,data$ncb),Pdc=matrix(rep(0.5,4),2,2),Pdn=rep(0.5,2)))
out<- mod$sample(
data=data,
init = inits,
iter_warmup = nb,
iter_sampling = ni,
chains=nc,
parallel_chains = nc,
refresh = 100)
out$save_object(file="Outputs/Koala_model_outputs.stan")
out$save_output_files(dir = "Outputs", basename = "KoalaData")
saveRDS(data, "Outputs/StanData.rds")
} else {
out<- readRDS("Outputs/Koala_model_outputs.stan")
}
Model Diagnostics
Based on the below MCMC chains, convergence and chain mixing occurred.
Show code
x <- out$summary(c("beta","Pdc","Pdn"))
write.csv(x, "Outputs/Parameter_estimates_final2.csv")
z <- out$summary(c("Nn_region","Np_region","Nhat_n","Nhat_p","Nhat"))
write.csv(z, file="Outputs/Abundance_estimates_final2.csv")
color_scheme_set()
mcmc_intervals(out$draws("beta")) + geom_vline(xintercept=0, linetype=2) +
theme_bw()
Show code
color_scheme_set("viridis")
mcmc_trace(out$draws("beta")) + theme_bw()
Model Predictions
Statewide Predictions
Below we map the koala densities across Victoria for forest classed as either native (> 5 ha native forest per square-kilometre) or plantation forest (> 1 ha plantation forest per square-kilometre). Areas are classed as plantation if the area of plantation is greater than the area of native forest for a 1km2 area.
Show code
# Only rerun if model is refreshed
## Native
preds<- out$summary(c("Nn"))
Nn <- out$draws(variable = "Nn", format = "draws_matrix")
trimed_mean_Nn <- apply(Nn, 2, function(x) {
mean(x, trim = 0.2)
})
preds$trimmedMean <- trimed_mean_Nn
saveRDS(preds, "Outputs/preds.rds")
native.rast <- raster(thresholdRaster$HabitatSuitability)
native.rast <- mask(native.rast, vic_bound)
native.rast[native_df$cell] <- preds$mean #native_df$NativeProp
# plot(native.rast)
writeRaster(native.rast, "Outputs/Abundance_prediction_native.tif", overwrite=TRUE)
native.se.rast <- native.rast
native.se.rast[native_df$cell] <- preds$sd/preds$mean
# plot(native.se.rast)
writeRaster(native.se.rast, "Outputs/SE_prediction_native.tif", overwrite=TRUE)
## Plantation
predsPlan<- out$summary(c("Np"))
Np <- out$draws(variable = "Np", format = "draws_matrix")
trimed_mean_Np <- apply(Np, 2, function(x) {
mean(x, trim = 0.2)
})
predsPlan$trimmedMean <- trimed_mean_Np
saveRDS(predsPlan, "Outputs/predsPlan.rds")
plantation.rast <- raster(thresholdRaster$HabitatSuitability)
plantation.rast <- mask(plantation.rast, vic_bound)
not.na <- which(is.na(values(plantation.rast))=='FALSE')
plantation.rast[not.na] <- 0
plantation.rast[plantation_df$cell] <- predsPlan$mean #plantation_df$PlantationProp
# plot(plantation.rast)
writeRaster(plantation.rast, "Outputs/Abundance_prediction_plantation.tif", overwrite=TRUE)
plantation.se.rast <- plantation.rast
plantation.se.rast[plantation_df$cell] <- predsPlan$sd/predsPlan$mean
# plot(plantation.se.rast)
writeRaster(plantation.se.rast, "Outputs/SE_prediction_plantation.tif", overwrite=TRUE)
Show code
# Discrete values prediction for native forest
preds<- readRDS("Outputs/preds.rds")
native.rast <- raster("Outputs/Abundance_prediction_native.tif")
coords<- as.data.frame(raster::coordinates(native.rast))
coords<- coords[native_df$cell,]
preds<- cbind(preds,coords[1:nrow(preds),], data.frame(NativeProp = native_df$NativeProp)) %>%
mutate(densitykm = mean*NativeProp)
preds$density<- preds$densitykm/100 #native_df$NativeProp[1:nrow(preds)]
#cp<- c(seq(0,3,0.5),4,20)
#labs<- c("0 - 0.5","0.5 - 1","1 - 1.5","1.5 - 2","2 - 2.5","2.5 - 3","3 - 4","4+")
cp<- c(seq(0,2,0.5),max(preds$density))
labs<- c("0 - 0.5","0.5 - 1","1 - 1.5","1.5 - 2","2+")
preds$density_disc<- cut(preds$density, cp, include.lowest=TRUE,right=FALSE)
p1<- vic_bound %>% ggplot() +
geom_sf(fill=NA) +
geom_tile(aes(x=x, y=y, fill=density_disc), data=preds) +
viridis::scale_fill_viridis(discrete=TRUE, direction=-1, option="D",labels=labs) +
geom_sf(fill=NA) +
labs(x="Longitude",y="Latitude",fill="Density (ha)", subtitle="A") +
theme_bw() +
theme(axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12),
plot.subtitle = element_text(face="bold", size=12),
legend.title = element_text(face="bold"))
# Discrete values prediction for plantations
predsPlan<- readRDS("Outputs/predsPlan.rds")
plantation.rast <- raster("Outputs/Abundance_prediction_plantation.tif")
coords<- as.data.frame(raster::coordinates(plantation.rast))
coords<- coords[plantation_df$cell,]
predsPlan<- cbind(predsPlan,coords[1:nrow(predsPlan),])
predsPlan$density<- predsPlan$mean #plantation_df$PlantationProp[1:nrow(predsPlan)]
cp<- c(seq(0,2,0.5), max(predsPlan$density))
labs<- c("0 - 0.5","0.5 - 1","1 - 1.5","1.5 - 2", "2+")
predsPlan$density_disc<- cut(predsPlan$density, cp, include.lowest=TRUE,right=FALSE)
p2<- vic_bound %>% ggplot() +
geom_sf(fill=NA) +
annotation_scale(location = "br", width_hint = 0.25) +
annotation_north_arrow(location = "br", which_north = "true", height=unit(0.5, "in"),
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
geom_tile(aes(x=x, y=y, fill=density_disc), data=predsPlan) +
scale_fill_viridis(discrete=TRUE, direction=-1, option="D") +
geom_sf(fill=NA) +
labs(x="Longitude",y="Latitude",fill="Density (ha)",subtitle="B") +
theme_bw() +
theme(axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12),
plot.subtitle = element_text(face="bold", size=12),
legend.title = element_text(face="bold"))
p1
Figure 6: Native forest koala densities
Show code
native.se.rast <- rast("Outputs/SE_prediction_native.tif") %>%
terra::mask(thresholdRaster$`KDE Distribution`)
native.se.rast.df <- as.data.frame(native.se.rast, xy = TRUE) %>%
sf::st_as_sf(., coords = c("x", "y"), crs = 3111) %>%
sf::st_join(delwp_region)
native.se.rast.df.avg <- native.se.rast.df %>%
sf::st_drop_geometry() %>%
group_by(DELWP_REGION) %>%
summarise(`CV mean` = mean(SE_prediction_native, na.rm = T),
`CV n` = n())
native.se.rast.df.g <- as.data.frame(native.se.rast, xy = TRUE) %>%
sf::st_as_sf(., coords = c("x", "y"), crs = 3111) %>%
sf::st_intersection(gips_bound %>% sf::st_transform(3111))
ggplot() +
geom_spatraster(data = native.se.rast, aes(fill = SE_prediction_native), na.rm = TRUE) +
geom_sf(data = vic_bound, fill = "grey", inherit.aes = FALSE, alpha = 0, colour = "grey70") +
# geom_sf(data = regions_3111, fill = "grey", inherit.aes = FALSE, alpha = 0.2, colour = "grey70") +
# geom_sf(data = statewide_water_simp, colour = "#2c7fb8", inherit.aes = FALSE, alpha = 0.5) +
# geom_tile(aes(fill = value), na.rm = TRUE) +
scale_fill_gradient2(na.value = "transparent", low = delwp_cols[["Environment"]], mid = delwp_cols[["Teal"]], high = delwp_cols[["Navy"]], trans = "log", breaks = c(0.5, 5, 50)) +
ylab("Latitude") +
xlab("Longitude") +
labs(fill = "CV") +
theme(axis.title = element_text(size = 40),
axis.text = element_text(size = 30),
legend.title = element_text(size=40),
legend.text = element_text(size=30),
legend.key.size = unit(6, 'cm')) +
theme_bw()
Figure 7: Coefficient of Variation for native forest koala densities
Regional Predictions
Based on the model we can also generate regional-specific predictions for koala abundance numbers in each DELWP region. These region-specific estimates are not directly comparable with the table in the previous report as a different habitat suitability threshold and native-area/plantation area classification is used.
Show code
NativeRegionalPred <- out$summary("Nn_region")
Nn_region <- out$draws(variable = "Nn_region", format = "draws_matrix")
trimed_mean_Nn_region <- apply(Nn_region, 2, function(x) {
mean(x, trim = 0.2)
})
NativeRegionalPred$trimmedMean <- trimed_mean_Nn_region
NativeRegionalPred_form <- NativeRegionalPred %>%
mutate(DELWP_REGION_CODE2 = 1:6) %>%
left_join(vic.regions.df %>%
dplyr::select(-DELWP_REGION_CODE), by = "DELWP_REGION_CODE2") %>%
dplyr::select(Region,
`Predicted Abundance` = mean,
`Standard Deviation` = sd,
`LCI [5%]` = q5,
`UCI [95%]` = q95) %>%
janitor::adorn_totals(., where = "row")
NativeRegionalPred_form %>%
kbl(digits = 0, format.args = list(big.mark = ","), caption = "Predicted Koala abundance in native forest and woodland across Victoria. The total state-wide prediction is provided, along with predictions for each of the six DELWP regions. LCI, lower credible interval; UCI, upper credible interval.") %>%
kable_styling(bootstrap_options = "striped") %>%
row_spec(row = 7, bold = TRUE)
| Region | Predicted Abundance | Standard Deviation | LCI [5%] | UCI [95%] |
|---|---|---|---|---|
| GIPPSLAND | 43,945 | 8,776 | 31,056 | 59,461 |
| PORT PHILLIP | 19,322 | 4,982 | 12,329 | 28,563 |
| HUME | 47,889 | 81,107 | 14,148 | 117,736 |
| BARWON SOUTH WEST | 199,901 | 38,449 | 141,084 | 266,766 |
| LODDON MALLEE | 6,820 | 2,653 | 3,615 | 11,399 |
| GRAMPIANS | 22,569 | 6,889 | 12,962 | 35,054 |
| Total | 340,446 | 142,855 | 215,193 | 518,980 |
East Gippsland Predictions
Based on the model we can obtain estimates specific to East Gippsland by cropping predictions to that area. Specifically, we define East Gippsland as the local government area (LGA) of East Gippsland.
Show code
Nn <- out$draws(variable = "Nn", format = "draws_matrix")
predEGsf <- preds %>%
sf::st_as_sf(., coords = c("x", "y"), crs = 3111) %>%
sf::st_intersection(gips_bound %>% sf::st_transform(3111))
ns <- predEGsf$variable
Nn_EG <- Nn[,ns]
Nn_EG_km <- as.matrix(Nn_EG) %*% diag(predEGsf[["NativeProp"]])
EGSum <- rowSums(Nn_EG_km)
PredMean <- apply(Nn_EG_km, 2, mean, trim = 0.2) # Should you trim some draws
# predEGsf_masked <- terra::extract(thresholdRaster$koalarecords_mask, vect(predEGsf))
# predEGsf_f <- predEGsf[which(predEGsf_masked$koalarecords_mask == 1),]
# predEG_f <- cbind(st_drop_geometry(predEGsf_f), as.data.frame(st_coordinates(predEGsf_f)))
predsEG <- bind_cols(st_drop_geometry(predEGsf), as.data.frame(st_coordinates(predEGsf)), data.frame(Trimmeddensitykm = PredMean))
egmedian <- mean(EGSum, trim = 0.2)
egLCI <- coda::HPDinterval(coda::as.mcmc(EGSum))[1,1]
egUCI <- coda::HPDinterval(coda::as.mcmc(EGSum))[1,2]
CVEG <- sd(EGSum)/mean(EGSum)
# f_egmedian <- sum(predEG_f$median)
# f_egvar <- sum(predEG_f$sd^2)
# f_egsd <- sqrt(f_egvar)
# f_egLCI <- max(c(f_egmedian - 1.96*f_egsd, 0))
# f_egUCI <- f_egmedian + 1.96*f_egsd
# cp<- c(seq(0,2,0.5),20)
# labs<- c("0 - 0.5","0.5 - 1","1 - 1.5","1.5 - 2","2+")
koala_map_proj <- cbind(st_drop_geometry(koala_map), as.data.frame(st_coordinates(koala_map %>% st_transform(3111)))) %>%
na.omit() %>%
bind_rows(cbind(data.frame(Presence = "Absent", as.data.frame(st_coordinates(GG_Gippsland_Survey_Data %>% st_intersection(gips_bound %>% st_transform(3111)))))))
gippspred_nocrop <- gips_bound %>%
sf::st_transform(3111) %>%
ggplot() +
geom_sf(fill=NA) +
geom_tile(aes(x=X, y=Y, fill=densityha), data=predsEG %>% mutate(densityha = densitykm/100)) +
viridis::scale_fill_viridis(discrete=FALSE, direction=-1, option="D") +
geom_sf(fill=NA) +
# geom_point(aes(x = X, y = Y, colour = Presence), shape = 19, size = 2, data = koala_map_proj) +
labs(x="Longitude",y="Latitude",fill="Density (ha)") +
theme_bw() +
theme(axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12),
plot.subtitle = element_text(face="bold", size=12),
legend.title = element_text(face="bold"))
# gippspred_crop <- gips_bound %>%
# sf::st_transform(3111) %>%
# ggplot() +
# geom_sf(fill=NA) +
# geom_tile(aes(x=X, y=Y, fill=density), data=predEG_f) +
# viridis::scale_fill_viridis(discrete=FALSE, direction=-1, option="D") +
# geom_sf(fill=NA) +
# geom_point(aes(x = X, y = Y, colour = Presence), shape = 19, size = 2, data = koala_map_proj) +
# labs(x="Longitude",y="Latitude",fill="Density (ha)") +
# theme_bw() +
# theme(axis.title.x = element_text(face="bold", size=12),
# axis.title.y = element_text(face="bold", size=12),
# plot.subtitle = element_text(face="bold", size=12),
# legend.title = element_text(face="bold"))
gippspred_nocrop
Figure 8: Model predictions for Koala density in East Gippsland. Model estimates include habitat according to a maxent threshold and withinthe smoothed density surface of the koala range, based on records in the previous 15 years.
Based on the Predictions from the point-process density model applied to a maxent threshold area the predicted trimmed mean density of koalas in East Gippsland is 4,537 [1,924,8,046].
If we compare the results to the last model the estimates for that were a median of 11,367 [10,981, 11,753].
Model Evaluation
Posterior Predictive Checks
We evaluate the models’ prediction power by conducting several posterior-predictive checks. These checks compare model posterior distributions to the observed value. If observed values fall within the central distribution of the models posterior distribution it is evidence that models have good predictive power.
Show code
ysc_rep <- out$draws(variable = c("ysc_rep"), format = "draws_matrix")
ydcd_rep<- out$draws(variable = c("ydcd_rep"), format = "draws_matrix")
ydcn_rep<- out$draws(variable = c("ydcn_rep"), format = "draws_matrix")
pred_rep<- cbind(ysc_rep,ydcd_rep,ydcn_rep)
obs<- c(data$ysc, rowSums(data$ydcd), rowSums(data$ydcn))
prop_zero<- function(x) mean(x == 0)
pp1<- ppc_stat(obs, pred_rep, stat=prop_zero, binwidth=0.005) +xlim(0.5, 0.75) +
ggtitle("Proportion zeros") + theme_bw()
pp2<- ppc_stat(obs, pred_rep, stat="mean") + ggtitle("Mean") + theme_bw()
pp3<- ppc_stat(obs, pred_rep, stat="sd") + ggtitle("Standard deviation") + theme_bw()
pp4<- ppc_stat(obs, pred_rep, stat="max") + xlim(0, 60) + ggtitle("Maximum") + theme_bw()
cowplot::plot_grid(pp1,pp2,pp3,pp4, labels = "AUTO")
Figure 9: Posterior predictive checks for the koala density model. Observed statistics that are compared to the posterior distributions are (A) proportion of zeros, (B) mean count, (C) standard deviation and (D) maximum count.
Additional figures show the posterior predictive checks at varying binned intervals of predicted values.
Show code
post1 <- pp_check(obs, pred_rep, fun="scatter_avg") + xlim(0, 30) + ylim(0, 30) + theme_bw()
post2 <- ppc_bars(obs, pred_rep, freq=FALSE) + theme_bw()
post3 <- ppc_rootogram(obs, pred_rep) + theme_bw()
post4 <- ppc_error_scatter_avg(obs, pred_rep)
cowplot::plot_grid(post1,post2,post3,post4, labels = "AUTO")
Figure 10: Additional posterior predictive scatter checks
We can also directly compare the koala counts obtained from the surveys to the predicted counts from the model for that particular survey. From this comparison we can see that the model predictions and survey observations are highly correlated.
Show code
par(mfrow=c(2,2),pty="s")
plot(data$ysc, colMeans(ysc_rep), xlim=c(0,30), ylim=c(0,30), ylab='Predicted', xlab='Observed',main="Single counts")
abline(a=0, b=1)
plot(rowSums(data$ydcd), colMeans(ydcd_rep), xlim=c(0,30), ylim=c(0,30), ylab='Predicted', xlab='Observed',
main="double counts - day")
abline(a=0, b=1)
plot(rowSums(data$ydcn), colMeans(ydcn_rep), xlim=c(0,30), ylim=c(0,30), ylab='Predicted', xlab='Observed',
main="double counts - night")
abline(a=0, b=1)
plot(obs, colMeans(pred_rep), xlim=c(0,30), ylim=c(0,30), ylab='Predicted', xlab='Observed', main="Combined")
abline(a=0, b=1)
Figure 11: Observed vs Predicted Koala counts across Victoria
Show code
med<- apply(pred_rep,2,median)
res<- DHARMa::createDHARMa(simulatedResponse=t(pred_rep),observedResponse = obs,
fittedPredictedResponse = med, integerResponse = TRUE)
plot(res)
Figure 12: QQ plots of model residuals
Show code
DHARMa::testUniformity(res)
Figure 13: Uniformity test of model residuals
One-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.086347, p-value = 1.068e-07
alternative hypothesis: two-sidedCovariate Effects
Show code
betas <- out$draws("beta", format = "draws_df")
names(betas) <- colnames(data$Xpn)
# head(betas)
# mean.betas <- colMeans(betas)
# betas_drawsdf <- as_draws(betas)
plot_stan <- function(betas, cov1, cov2, background) {
mean.betas <- colMeans(betas)
mincov1 <- min(background_df[[cov1]], na.rm = TRUE)
maxcov1 <- max(background_df[[cov1]], na.rm = TRUE)
cov1_seq <- seq(mincov1, maxcov1, 0.1)
mincov2 <- quantile(background_df[[cov2]], na.rm = TRUE ,0.05)
maxcov2 <- quantile(background_df[[cov2]], na.rm = TRUE, 0.95)
cov2_seq <- seq(mincov2, maxcov2, 0.1)
cov.pred <- matrix(nrow=length(cov1_seq), ncol=length(cov2_seq))
comb_var <- which(stringr::str_detect(colnames(betas), paste(cov1,cov2, sep=":")))
for(i in 1:nrow(cov.pred)){
for(j in 1:ncol(cov.pred)){
cov.pred[i, j] <- exp(mean.betas[["(Intercept)"]] + mean.betas[[cov1]]*cov1_seq[i] + mean.betas[[cov2]]*cov2_seq[j] +
mean.betas[[comb_var]] * (cov1_seq[i] * cov2_seq[j]))
}
}
cov.pred.df <- as.data.frame(cov.pred) %>%
`colnames<-`(as.character(cov2_seq))
cov.pred.df[[cov1]] <- cov1_seq
cov.pred.df.melt <- as.data.table(cov.pred.df) %>%
melt.data.table(id.vars = cov1, variable.name = cov2) %>%
dplyr::mutate(value = case_when(value > 5 ~ 5,
TRUE ~ value))
cov.pred.df.melt[[cov2]] <- as.numeric(as.character(cov.pred.df.melt[[cov2]]))
value <- "value"
ggplot(cov.pred.df.melt) +
geom_tile(aes_string(x = cov1, y = cov2, fill = value)) +
theme_minimal() +
scale_fill_gradientn(colours = c(scales::muted("#b2182b"),"#f7f7f7","#2166ac")) +
labs(fill = "Impact on \nDensity")
# fields::image.plot(cov1_seq, cov2_seq, cov.pred, axis.args=list(title=c(cov1, cov2)))
}
cowplot::plot_grid(
plot_stan(betas, "Rainfall", "HighPrefEucNBR", background = background_df),
plot_stan(betas, "Rainfall", "MedPrefEucNBR", background = background_df),
plot_stan(betas, "Rainfall", "LowPrefEucNBR", background = background_df),
plot_stan(betas, "Rainfall", "PlantationEucNBR", background = background_df),
plot_stan(betas, "Altitude", "HighPrefEucNBR", background = background_df),
plot_stan(betas, "Altitude", "MedPrefEucNBR", background = background_df),
plot_stan(betas, "Altitude", "LowPrefEucNBR", background = background_df),
plot_stan(betas, "Altitude", "PlantationEucNBR", background = background_df),
plot_stan(betas, "MTWQ", "HighPrefEucNBR", background = background_df),
plot_stan(betas, "MTWQ", "MedPrefEucNBR", background = background_df),
plot_stan(betas, "MTWQ", "LowPrefEucNBR", background = background_df),
plot_stan(betas, "MTWQ", "PlantationEucNBR", background = background_df), ncol = 3)
Figure 14: Impact of model covariates on koala density. Interaction effects of habitat covariates on the predicted density of koalas. Estimates are limited/truncated at 5 to minimises extrapolation into covariate spaces not present within the sampled data.
Show code
quants <- c(0.50,0.05,0.95)
BetaSummary <- apply(betas[,1:22] , 2 , quantile , probs = quants , na.rm = TRUE )
Beta_df <- as.data.frame(BetaSummary) %>% t()
Beta_df %>%
kbl(caption = "Parameter estimates for the top Koala abundance model, incorporating imperfect detection. Only estimates of the fixed effects are shown. Variables were standardised prior to model fitting, meaning parameter estimates are directly comparable. Feed tree preference layers are the average neighbourhood prevalence in all cases.", digits = 2) %>%
kable_styling(bootstrap_options = "striped")
| 50% | 5% | 95% | |
|---|---|---|---|
| (Intercept) | -2.48 | -3.35 | -1.58 |
| hab.valsNative | 0.02 | -0.27 | 0.30 |
| Rainfall | -0.35 | -0.82 | 0.11 |
| Altitude | 0.04 | -0.36 | 0.44 |
| MTWQ | -1.34 | -2.09 | -0.67 |
| Slope | -0.12 | -0.40 | 0.14 |
| HighPrefEucNBR | -0.94 | -1.73 | -0.17 |
| MedPrefEucNBR | 0.04 | -0.25 | 0.35 |
| LowPrefEucNBR | -0.63 | -1.62 | 0.24 |
| PlantationEucNBR | -0.78 | -1.10 | -0.48 |
| Rainfall:HighPrefEucNBR | 0.29 | -0.10 | 0.66 |
| Rainfall:MedPrefEucNBR | 0.96 | 0.58 | 1.36 |
| Rainfall:LowPrefEucNBR | -0.41 | -0.86 | 0.03 |
| Rainfall:PlantationEucNBR | -0.20 | -0.33 | -0.07 |
| Altitude:HighPrefEucNBR | -0.63 | -0.98 | -0.30 |
| Altitude:MedPrefEucNBR | 0.00 | -0.30 | 0.30 |
| Altitude:LowPrefEucNBR | -0.10 | -0.53 | 0.34 |
| Altitude:PlantationEucNBR | -0.33 | -0.45 | -0.21 |
| MTWQ:HighPrefEucNBR | -0.83 | -1.50 | -0.18 |
| MTWQ:MedPrefEucNBR | 0.72 | 0.37 | 1.08 |
| MTWQ:LowPrefEucNBR | 0.03 | -0.89 | 0.88 |
| MTWQ:PlantationEucNBR | -0.76 | -1.08 | -0.45 |
Detection Probabilities
Show code
det_prob <- out$summary(c("Pdc","Pdn"))
det_prob$Time <- c(rep("Diurnal", 4), rep("Nocturnal", 2))
det_prob$Habitat <- c("Plantation", "Native", "Plantation", "Native", "Native", "Native")
det_prob$Observer <- c(1,1,2,2,1,2)
det_prob %>%
arrange(Time, Habitat) %>%
dplyr::select(Time, Habitat, Observer, `Mean Det` = mean, SD = sd, `5 % LCI` = q5, `95 % UCI` = q95) %>%
kbl(digits = 2) %>%
kable_styling(bootstrap_options = "striped") %>%
collapse_rows(1:3)
| Time | Habitat | Observer | Mean Det | SD | 5 % LCI | 95 % UCI |
|---|---|---|---|---|---|---|
| Diurnal | Native | 1 | 0.80 | 0.02 | 0.77 | 0.83 |
| 2 | 0.69 | 0.02 | 0.66 | 0.73 | ||
| Plantation | 1 | 0.68 | 0.02 | 0.63 | 0.71 | |
| 2 | 0.67 | 0.02 | 0.63 | 0.71 | ||
| Nocturnal | Native | 1 | 0.19 | 0.13 | 0.04 | 0.45 |
| 2 | 0.08 | 0.07 | 0.01 | 0.22 |
Fire Impacts
Based on the predictions in East Gippsland we can overlay the fire severity map to visualize the number of koalas impacted by fire. Largely areas of predicted koala density occur outside high severity fire areas.
Show code
f_sev <- extract(fire_severity, vect(predEGsf))
names(f_sev) <- c("cell", "Fire Severity")
# table(f_sev$`Fire Severity`)
predEGsf_fire_sev <- bind_cols(predsEG, f_sev)
predEGsf_fire_sev_class <- predEGsf_fire_sev %>%
group_by(`Fire Severity`) %>%
summarise(`Koalas Impacted` = round(sum(densitykm)))
gips_bound %>%
sf::st_transform(3111) %>%
ggplot() +
geom_sf(fill=NA) +
tidyterra::geom_spatraster(data = fire_severity_c, na.rm = TRUE, inherit.aes = FALSE, alpha = 0.8) +
geom_tile(aes(x=X, y=Y), data=predsEG %>% filter(mean > 0), fill="Blue", alpha = 0.75) +
scale_fill_gradientn(
colors = hcl.colors(10, "Reds", rev = TRUE),
na.value = NA, guide = guide_colourbar(title = "Fire \nSeverity\n")) +
geom_sf(fill=NA) +
# geom_point(aes(x = X, y = Y, colour = Presence), shape = 19, size = 2, data = koala_map_proj) +
labs(x="Longitude",y="Latitude",fill="Density (ha)") +
theme_bw() +
theme(axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12),
plot.subtitle = element_text(face="bold", size=12),
legend.title = element_text(face="bold"))
Figure 15: Fires severity and areas where predicted koala density is non-zero (trimmed mean > 0)
Show code
predEGsf_fire_sev_class %>%
kbl(caption = "Predicted koala locations in East Gippsland largely fall within areas with little or no fire severity from the 2020 bushfires") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE)
| Fire Severity | Koalas Impacted |
|---|---|
| 0 | 1919 |
| 1 | 184 |
| 2 | 225 |
| 3 | 532 |
| 4 | 111 |
| 5 | 823 |
| 6 | 447 |
| NA | 485 |
Limitations and Future Directions
By updating the previous model (with several minor changes) and new double count data from East Gippsland (165 sites from Greater Glider and Koala Surveys) we were able to compare previous distributions and densities to regenerated ones. However given the large uncertainty in predictions, especially for East Gippsland, there are several limitations in the predictions made here:
Due to the predominant amount of koala detections being in koala hotpots (e.g. South-West region and Otways), similar habitats in Eastern Victoria may be modelled as having high koala densities despite no previous records nearby. To potentially improve these spatial variations, future models may be better placed to include spatial correlation estimates between sites (e.g. Spatial Gaussian processes).
More balanced surveying across Victoria will help counteract regional bias in the data.
More distance-sampling to feasibly have a large enough sample to begin include detection functions within the analysis.
Consider the use of two-step occupancy and abundance models if there is evidence that koala occupancy may be a different process than koala abundance.
Predictions have only been made to area that has been included in a maxent threshold, and contains a given proportion of native vegetation. These restrictions can have a large impact on the total number of koalas predicted across regions and the state.
Show code
dshare <- here::here("_posts", "2022-04-06-koala-abundance-monitoring", "data", "DataShare")
sf::write_sf(koala_transects, paste0(dshare, "/koala_transects.kml"))
readr::write_csv(raw_headers, paste0(dshare, "/raw_headers.csv"))
readr::write_csv(raw_animals, paste0(dshare, "/raw_animals.csv"))
readr::write_csv(all_transects, paste0(dshare, "/all_transects.csv"))
readr::write_csv(all_koala_observations, paste0(dshare, "/all_koala_observations.csv"))
Aiello-Lammens, Matthew E., Robert A. Boria, Aleksandar Radosavljevic, and Bruno Vilela. 2015. “spThin: An R Package for Spatial Thinning of Species Occurrence Records for Use in Ecological Niche Models.” Ecography 38 (5): 541–45. https://onlinelibrary.wiley.com/doi/10.1111/ecog.01132.
Department of Environment, Water & Planning, Land. 2022. “Fire Severity Map of the Major Fires in Gippsland and North East Victoria in 2019/20 (Version 1.0).” 2022. http://services.land.vic.gov.au/catalogue/metadata?anzlicId=ANZVI0803008638&publicId=guest&extractionProviderId=1.
Stevenson, Matilda, Martin Westgate, and Peggy Newman. 2022. Galah: Atlas of Living Australia (ALA) Data and Resources in r. https://CRAN.R-project.org/package=galah.
Vignali, Sergio, Arnaud G. Barras, Raphaël Arlettaz, and Veronika Braunisch. 2020. “SDMtune: An r Package to Tune and Evaluate Species Distribution Models.” Ecology and Evolution 10: 11488–506. https://doi.org/10.1002/ece3.6786.