07-projections.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

<!-- 
This file is part of a gitbook that should be cited as: 

Valle, M., Citores, L., Ibaibarriaga, L., Chust, C. (2023) GAM-NICHE: Shape-Constrained GAMs to build Species Distribution Models under the ecological niche theory. AZTI. https://doi.org/10.57762/fzpy-6w51 

This tutorial has been supported by the European Union’s Horizon 2020 research and innovation programme under grant agreements No 862428 MISSION ATLANTIC project
-->

# Prediction and maps

In this chapter we predict from the fitted model and produce final SDMs maps.

First, we load a list of required libraries.
```{r, eval=T,message=FALSE,warning=FALSE}
requiredPackages <- c(
  "here", 
  "rstudioapi", 
  "ggplot2", 
  "tidyverse", 
  # "rgdal", 
  "raster", 
  "maps", 
  "RColorBrewer", 
  "scam", 
  "ggpubr"
  )
```

We run a function to install the required packages that are not in our system and load all the required packages.
```{r, eval=T,message=FALSE,warning=FALSE}
install_load_function <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

install_load_function(requiredPackages)

```

We define some overall settings.
```{r, eval=T,message=FALSE,warning=FALSE}
# General settings for ggplot (black-white background, larger base_size)
theme_set(theme_bw(base_size = 16))
```

## Prepare environmental data

In previous steps (see Chapter 2), we have defined the study area that defines the extent of our spatial data. We load the `study_area` object that is a SpatialPolygonsDataFrame class:
```{r, eval=T, warning=F, message=F}
load(here::here ("data", "spatial", "study_area.RData"))
```

And we load the rasterStack with the downloaded environmental data.
```{r, eval=T, warning=F, message=F}
mylayers<-stack(here::here ("data", "env", "mylayers.tif"))
```

We transform the environmental data set first into a data frame, and then into a SpatialDataFrame.
```{r, eval=T, warning=F, message=F}
env_dataframe <- raster::as.data.frame(mylayers, xy=TRUE)

summary(env_dataframe)

names(env_dataframe) <- c("x", "y", "BO2_chlomean_ss", "BO2_salinitymean_ss", "BO_damean", "BO_sstmean")
```

## Projection

We load the selected model and predict into the whole environmental data. 
```{r, eval=T, warning=F, message=F}
# Load SC-GAM model
load(here::here("models", "selected_model.Rdata"))

# Predicting 
predict <- predict(selected_model,newdata=env_dataframe,type ="response",se.fit=T)         

env_dataframe$fit<-predict$fit
env_dataframe$se.fit<-predict$se.fit

save(env_dataframe, file="results/projection.Rdata")
```

## Mapping

```{r, eval=T, warning=F, message=F}
# Load PA data
load(here::here ("data", "outputs_for_modelling", "PAdata_with_env.Rdata"))


proj_map <-ggplot()+
  geom_raster(data=subset(env_dataframe),
              aes(x,y,fill=fit)) +
  scale_fill_gradient2(low="blue", 
                       mid="orange",
                       high="red",
                       midpoint = 0.5,
                       limits = c(0,1)) +
  ggtitle("Occurrence probabilty Thunnus alalunga")+ 
  geom_point(data=subset(data,occurrenceStatus==1),
             aes(LON,LAT),
             col=1,
             size=0.3) +
  theme_pubclean(base_size = 14)+
  theme(panel.background = element_blank(),
        plot.title = element_text(face = "italic"), 
        #text = element_text(size = 14), 
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        legend.position="right") +
  labs(y="latitude", x = "longitude")
  
print(proj_map)
```

We finally save the projection map.
```{r, eval=T, warning=F, message=F}
ggsave(filename= "Thunnus_alalunga_proj_map.tif", 
        plot=proj_map, 
        device="tiff",
        path=here::here ("plots", "projections"), 
        height=22, width=30,
        units="cm", dpi=300)
```