02-data_pa.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
-->

# Presence-absence data

In this chapter we first, download occurrence data from global open-access datasets such as Global Biodiversity Information Facility (GBIF, https://www.gbif.org/) and Ocean Biodiversity Information System (OBIS, https://obis.org/); second, clean downloaded data reformating, renaming fields and removing outliers data; and lastly, we generate a set of pseudoabsence points along the defined study area. 

First we load a list of required libraries.
```{r, eval=T,message=FALSE,warning=FALSE}
requiredPackages <- c(
  "here", 
  "rstudioapi",  
  "ggplot2", 
  "robis", 
  "rgbif",
  "CoordinateCleaner", 
  "sf", 
  "data.table", 
  "dplyr", 
  "tidyr", 
  "marmap", 
  "tidyverse", 
  "scales",
  "ggridges",
  "maps",     
  "mapdata",
  "mapproj",
  "mapplots",
  "gridExtra",
  "lubridate",
  "raster" 
    )
```

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))
```

## Download presence data

In this section we download presence data from global public datasets. 

To do so, we first define a study area, in this case we select the Atlantic Ocean based on the The Food and Agriculture Organization (FAO) Major Fishing Areas for Statistical Purposes and we remove Black Sea subarea.   
```{r, eval=T,message=FALSE,warning=FALSE}
# we could download the shapefile with the FAO fishing areas from the url where FAO shapefile is stored
# uncommenting the next two lines: 

# url<-"https://www.fao.org/fishery/geoserver/area/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=area%3AFAO_AREAS&maxFeatures=50&outputFormat=SHAPE-ZIP"
# download.file(url,"data/spatial/FAO_AREAS.zip",mode="wb")

# Unzip downloaded file
unzip(here::here ("data", "spatial", "FAO_AREAS.zip"),exdir="data/spatial")

# Load FAO (spatial multipolygon)
FAO<- st_read(file.path("data", "spatial", "FAO_AREAS.shp"))

# Select Atlantic Ocean FAO Area 
FAO_Atl <- FAO[FAO$OCEAN=="Atlantic",]

# Select Black Sea subarea 
Black_Sea <- FAO_Atl[FAO_Atl$ID=="20",]

# Transform to sf objects
FAO_Atl.sf <- st_as_sf(FAO_Atl)
Black_Sea.sf <- st_as_sf(Black_Sea)

# Remove Black sea using st_difference (reverse of st_intersection)
FAO_Atl_no_black_sea <- st_difference(FAO_Atl.sf,Black_Sea.sf) %>%   dplyr::select (F_AREA)

# Transform to spatial polygons dataframe
study_area<- sf:::as_Spatial(FAO_Atl_no_black_sea)

plot(study_area)

# Remove unused files
rm(FAO, FAO_Atl, FAO_Atl.sf, FAO_Atl_no_black_sea, Black_Sea, Black_Sea.sf)
```

Download occurrence data from OBIS and GBIF using scientific name. 

In this case we select Albacore tuna species (<em>Thunnus alalunga</em>). This can take a long time, so we upload data downloaded previously.

```{r, eval=T, message=FALSE, warning=FALSE}

# To get data from OBIS
# mydata.obis<-robis::occurrence(scientificname="Thunnus alalunga")

# To get data from GBIF
# mydata.gbif<-occ_data(scientificName="Thunnus alalunga", hasCoordinate = TRUE, limit=100000)$data

# Given that this can take a long time, we upload data downloaded previously
load(here::here ("data", "occurrences", "mydata.obis.RData"))
load(here::here ("data", "occurrences", "mydata.gbif.RData"))
```


We now check the downloaded data and select the fields of interest.
```{r, eval=T,message=FALSE,warning=FALSE}
# Check names for GBIF data
names(mydata.gbif)
  
# Select columns of interest
mydata.gbif <- mydata.gbif %>%
                dplyr::select("acceptedScientificName",
                  "decimalLongitude",
                  "decimalLatitude",
                  "year",
                  "month",
                  "day",
                  "eventDate",
                  "depth")

# Check names in for OBIS data
names(mydata.obis)

# Select columns of interest
mydata.obis <-  mydata.obis %>%
                dplyr::select("scientificName",
                  "decimalLongitude",
                  "decimalLatitude",
                  "date_year",
                  "month",
                  "day",
                  "eventDate",
                  "depth",
                  "bathymetry",
                  "occurrenceStatus",
                  "sst")

```

Reformat the data adding a new field and renaming some columns from mydata.gbif dataframe in order to have the same columns and be able to join both downloaded datasets.
```{r, eval=T,message=FALSE,warning=FALSE}
mydata.gbif <- mydata.gbif %>% 
    dplyr::rename(scientificName= "acceptedScientificName") %>% 
    dplyr::rename(date_year = "year") %>% 
    dplyr::mutate(bathymetry= NA) %>% 
    dplyr::mutate(occurrenceStatus=1) %>% 
    dplyr::mutate(sst= NA)

# Join data from OBIS and GBIF 
mydata.fus<-rbind(mydata.obis,mydata.gbif)
  
# Assign unique scientific name 
mydata.fus <- mydata.fus %>% 
    dplyr::mutate(scientificName= paste(mydata.obis$scientificName[1]))

# Remove unused files
rm(mydata.gbif, mydata.obis)
```

We now clean downloaded raw data.
```{r, eval=T,message=FALSE,warning=FALSE}
# Give date format to eventDate and fill out month and date_year columns
  mydata.fus$eventDate <- as.Date(mydata.fus$eventDate)
  mydata.fus$date_year <- as.numeric(mydata.fus$date_year)
  mydata.fus$month <- as.numeric(mydata.fus$month)

# Assign 1 value to occurrenceStatus
mydata.fus <- mydata.fus  %>% 
    dplyr::mutate(occurrenceStatus = 1)
```

We remove outliers based on distance method (total distance= 1000 km) available in the `CoordinateCleaner` package.
```{r, eval=T,message=FALSE,warning=FALSE}

out.dist <- cc_outl(x=mydata.fus,
                lon = "decimalLongitude", lat = "decimalLatitude",
                species = "scientificName",
                method="distance", tdi=1000,
                thinning=T, thinning_res=0.5,
                value="flagged") 

# Remove outliers from the data
mydata.fus <- mydata.fus[out.dist, ]
```

Remove duplicates.
```{r, eval=T,message=FALSE,warning=FALSE}
# First create a vector containing longitude, latitude and event date information
date <- cbind(mydata.fus$decimalLongitude,mydata.fus$decimalLatitude,mydata.fus$eventDate)

# Remove the duplicated records
mydata.fus<-mydata.fus[!duplicated(date),]

# Remove unused files
rm(date)
```

Mask retrieved occurrence data to our study area and add bathymetry value to each occurrence point.
```{r, eval=T,message=FALSE,warning=FALSE}
# Assign coordinate format and projection to be able to use FAO Atlantic as a mask
dat <- data.frame(cbind(mydata.fus$decimalLongitude,mydata.fus$decimalLatitude))
ptos<-as.data.table(dat,keep.columnnames=TRUE)
  
coordinates(ptos) <- ~ X1 + X2

# Assign projection
proj4string(ptos) <-proj4string(study_area)
  
# Select only occurrences from FAO Atlantic
match2<-data.frame(subset(mydata.fus,!is.na(over(ptos, study_area)[,1])))
  
# Extract the FAO area of each point
match3<-data.frame(subset(over(ptos, study_area), !is.na(over(ptos, study_area)[,1])))
  
# Create data frame with area, name, long, lat and year 
df0<-cbind(F_AREA=match3$F_AREA,match2)[,c("F_AREA","scientificName","decimalLongitude","decimalLatitude","date_year","occurrenceStatus")]

# Rename some columns
names(df0)[3:5]<-c("LON","LAT","YEAR")

options(timeout=600)

# Add bathymetry from NOAA
bathy <- marmap::getNOAA.bathy(lon1=-100,lon2=30,lat1=-41,lat2=55, resolution = 1,keep=TRUE, antimeridian=FALSE, path= here::here("data", "spatial"))

# instead, we could upload a file already downloaded
# load(here::here("data", "spatial","bathy.RData")) 

df0$bathymetry <- get.depth(bathy, df0[,c("LON","LAT")], locator=F)$depth

# Remove unused files
rm(mydata.fus, match2, match3, dat, ptos)
```

We plot the occurrence data.
```{r, eval=T,message=FALSE,warning=FALSE}
ggplot() +
   geom_path(data = study_area, 
             aes(x = long, y = lat, group = group),
             color = 'gray', linewidth = .2) +
   geom_point(data=df0, aes(x=LON, y=LAT)) 

```

And, we save the data.
```{r, eval=T,message=FALSE,warning=FALSE}
save(df0, file = "data/occurrences/occ.RData")
save(study_area, file = "data/spatial/study_area.RData")
```

## Create pseudo-absence data

After saving presence data for the species of interest we need to generate absence information in order to work with logistic regression SDM afterwards. In this case we create pseudo-absence data with a constant prevalence of 50% [@mcpherson_2004]; [@barbetmassin_etal_2012].

For that, we generate a buffer around each presence data point, with a radius of 100km, where no points can be generated. 

Before starting the pseudo-absence generation process, we delete observations in land (positive bathymetry) and select data between 2000 and 2014 (same temporal range as the environmental data that we will download later).
```{r,  eval=T,message=FALSE,warning=FALSE}
# Remove points in land
df0<-subset(df0,bathymetry<0)

# Select only years from 2000 to 2014
df0<- subset(df0, YEAR<=2014 & YEAR>=2000)

```

Then we transform the presence data frame to "SpatialPointsDataFrame" class object and to "sf" class object so we can operate and plot easily with spatial data tools. We can see the characteristics of the object through its summary.
```{r,  eval=T,message=FALSE,warning=FALSE}
# Convert to spatial point data frame
df<-df0 ; coordinates(df)<- ~LON+LAT
crs(df)<-crs(study_area)

# Convert to sf
df.sf<-st_as_sf(df)
study_area.sf<-st_union(st_as_sf(study_area))


summary(df)
```

We can easily plot sf objects using ggplot, here we plot the area of study and the presence data points.
```{r,  eval=T,message=FALSE,warning=FALSE}
ggplot(study_area.sf) + 
  geom_sf() + 
  geom_sf(data=st_union(df.sf),
          size=1,
          alpha=0.5)
```

We save a base plot (p0) with the world map and our area of study in blue.
```{r, eval=T,message=FALSE,warning=FALSE}
# Basic ggplot
global <- map_data("worldHires")

p0 <- ggplot() + 
  annotation_map(map=global, fill="grey")+
  geom_sf(data=study_area.sf,fill=5)

print(p0)
```

In order to generate a buffer around each occurrence data point, we need to work with euclidean distances, so first, we need to transform the decimal latitude and longitude values to UTM.
```{r, eval=T,message=FALSE,warning=FALSE}
# Function to find your UTM. 
lonlat2UTM = function(lonlat) {
  utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1
  if(lonlat[2] > 0) {
    utm + 32600
  } else{
    utm + 32700
  }
}

(EPSG_2_UTM <- lonlat2UTM(c(mean(df$LON), mean(df$LAT))))

# Transform study_area and data points to UTMs (in m)
aux <- st_transform(study_area.sf, EPSG_2_UTM)
df.sf.utm <- st_transform(df.sf, EPSG_2_UTM)
```

Now, we can create buffers of 100 km around the points and join the resulting polygons. Then this buffer is intersected with the area of study defining the area where the pseudo-absences can be generated. To visualize the defined areas, we plot the buffers in red and the area that we will use to generate pseudo-absences in green.
```{r, eval=T,message=FALSE,warning=FALSE}
# Create buffers of 100000m
buffer <- st_buffer(df.sf.utm, dist=100000)
buffer <- st_union(buffer)

# Intersect the are with the buffer
aux0 <- st_difference(aux, buffer)

# ggplot for all data
p0 +
  geom_sf(data=aux0,fill=3) +
  geom_sf(data=buffer,fill=2)

#zoom
p0 +
  geom_sf(data=aux0,fill=3) +
  geom_sf(data=buffer,fill=2)+
  coord_sf(xlim=c(-95,-82), ylim=c(22,31))
```

We create a data frame for pseudo-absences with the same dimensions as the presences data frame.
```{r, eval=T,message=FALSE,warning=FALSE}
# Generate the pseudo-absence data frame
pseudo <- matrix(data=NA, nrow=dim(df0)[1], ncol=dim(df0)[2])
pseudo <- data.frame(pseudo)
names(pseudo) <- names(df0)
```

To generate the pseudo-absence data points, we sample randomly from the defined area and we extract their latitude and longitude to incorporate them in the created data frame. We set the occurrenceStatus equal to 0 as they are absences.
```{r, eval=T,message=FALSE,warning=FALSE}
# Set the seed
set.seed(1)

# Sample from the defined area
rp.sf <- st_sample(aux0, size=dim(df.sf.utm)[1], type="random") # randomly sample points

# Transform to lat and lon and extract coordinates as data.frame
rp.sf <- st_transform(rp.sf, 4326)
rp <- as.data.frame(st_coordinates(rp.sf)) 
pseudo$LON <- rp$X
pseudo$LAT <- rp$Y

# Complete the rest of columns
pseudo$scientificName <- df0$scientificName
pseudo$occurrenceStatus  <- 0
```

We can plot the generated pseudo-absence data (in pink) in the map, together with the presence data points (in black).
```{r, eval=T,message=FALSE,warning=FALSE}
p0 +
  geom_sf(data=rp.sf, col=6, shape=4,size=0.5)+
  geom_sf(data=df.sf.utm, col=1, alpha=0.8,size=0.5)+
  ggtitle(unique(df$scientificName))

# Zoom
p0 +
  geom_sf(data=rp.sf, col=6, shape=4,size=1)+
  geom_sf(data=df.sf.utm, col=1, alpha=0.8,size=0.5)+
  coord_sf(xlim=c(-95,-82), ylim=c(23,31))+
  ggtitle(unique(df$scientificName))
```

Finally we join the presence and pseudo-absence data frames selecting the columns of interest and save the new data frame.
```{r, eval=T,message=FALSE,warning=FALSE}
# Join the two data sets 
PAdata <- rbind(df0, pseudo)[,c("scientificName","LON","LAT","YEAR","occurrenceStatus")]

# Save the final dataset of occurrence and pseudo-absence points
save(list=c("PAdata"),file=file.path("data","outputs_for_modelling",file="PAdata.RData"))
```