# Example of usage of the Chesapeake Bay atlas (https://doi.org/10.17882/99441)
# MACAN webinar of 2024-09-05, pst-laurent@vims.edu
#
# A map of habitat suitability for sandbark shark (Carcharhinus plumbeus) is
#   generated for the month of August.
 
# Import the different libraries that we will need. If a library/package isn't
#   already installed on your computer, type:
#   install.packages( 'name_of_library' )

  library( ncdf4        ) # To read NetCDF files.
  library( lattice      ) # To generate figures (levelplot).
  library( RColorBrewer ) # To access a collection of color palettes.
 
# Read from the atlas the different variables that we will need. We specifically
#   select bottom fields since this is where the shark lives.

  atla <- nc_open( 'atlas_chesbay_v20240319.nc' )

  long <- ncvar_get( atla, 'longitude'          ) # Dims: 348
  lati <- ncvar_get( atla, 'latitude'           ) # Dims:     567
  salt <- ncvar_get( atla, 'salinity_bottom'    ) # Dims: 348x567x12 (PSU    ).
  temp <- ncvar_get( atla, 'temperature_bottom' ) # Dims: 348x567x12 (deg.C  ).
  dio2 <- ncvar_get( atla, 'O2_bottom'          ) # Dims: 348x567x12 (mg-O2/L).

  nc_close( atla )                                # Close file once done.
  rm(       atla )                                # No longer necessary.

# Read the Generalized Additive Models (GAMs) from Crear et al. 2020
#   ( https://doi.org/10.3354/meps13483 ) for dissolved O2, practical salinity,
#   and potential temperature (their Figure 1). First column is environmental
#   variable, second column is the associated metric of suitability.

  gao2 <- read.csv( 'gam_crear2020_o2.csv'   ) # Dims: 43x2
  ga_s <- read.csv( 'gam_crear2020_salt.csv' ) # Dims: 24x2
  ga_t <- read.csv( 'gam_crear2020_temp.csv' ) # Dims: 21x2

# The original GAMs are restricted to a certain range of environmental
#   conditions. Extend the ranges with plausible values to cover all possible
#   conditions:

  gao2 <- rbind( c(0, 0.22), gao2,     c(12, 0), c(20, 0)           )
  ga_s <- rbind( c(0, 0.00), c( 3, 0), ga_s,     c(40, 0), c(50, 0) )
  ga_t <- rbind( c(0, 0.00), c(10, 0), ga_t,     c(35, 0)           )

# Convert the environmental conditions into metrics of habitat Suitability using
#   linear interpolation:

  s_o2 <- approx( gao2[,1], gao2[,2], dio2 ) # Suitability wrt O2.
  s__s <- approx( ga_s[,1], ga_s[,2], salt ) # Suitability wrt Salinity.
  s__t <- approx( ga_t[,1], ga_t[,2], temp ) # Suitability wrt Temperature.
  rm( gao2, ga_s, ga_t ) # No longer necessary.

# By default, the command 'approx' generates a 'data frame' containing both x,y
#   values. Extract the y values and reshape them into an array of same
#   dimensions as the environmental conditions (348x567x12):

  s_o2 <- array( data = s_o2$y, dim = dim( dio2 ) )
  s__s <- array( data = s__s$y, dim = dim( salt ) )
  s__t <- array( data = s__t$y, dim = dim( temp ) )

# Habitat Suitability Index taking the 3 environmental conditions into account
#   simultaneously:

  hsi  <- s_o2 * s__s * s__t
  rm( s_o2, s__s, s__t ) # No longer necessary.

# Select the MONTh of August, specifically.

  mont <- 8 # August;

# Generate a grid in the format expected by function 'levelplot':

  grid <- expand.grid( lon = long, lat = lati )

# Define the Contour INTervals for our figure of HSI: Zero to one, by intervals
#   of 0.01:

  cint = seq( from = 0, to = 1, by = 0.01 )

# Define a ColorMAP using the 'Spectral' palette; its length should be one less
#   than the length of the 'cint' vector (like the gaps between the teeth of a
#   fork):

  cmap <- colorRampPalette( brewer.pal( 11, 'Spectral') )( length( cint ) - 1 )

# Define the geographical extent of the figure. The values below correspond
#   approximately to the full atlas but one can focus on a specific sub-region
#   of the Bay by adjusting these values:

  xlim <- c( - 77.4, - 75.1 ) # Longitudes (full extent of atlas)
  ylim <- c(   36.6,   39.6 ) # Latitudes  (full extent of atlas)
# xlim <- c( - 76.5, - 76.1 ) # Longitudes (Elizabeth River     )
# ylim <- c(   36.7,   37.2 ) # Latitudes  (Elizabeth River     )

# Generate the HSI figure using 'levelplot'. Recall that hsi has dimensions
#   348x567x12 with the third dimension corresponding to the MONth of the year.
#   Argument 'at' informs levelplot to use the Contour INTervals defined above.
#   Argument 'col.regions' provide to levelplot our ColorMAP defined above.
#   Arguments 'xlim','xlim' define the geographical extent of the figure.

  levelplot( hsi[,,mont] ~ lon * lat, data = grid, at = cint, pretty = T, col.regions = rev( cmap ), xlim = xlim, ylim = ylim )