# 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 the calcification rate of Eastern oyster (Crassostrea virginica) is
#   generated for two contrasting months of the year, January and September.
 
# Import the different libraries that we will need. If they aren't already
#   installed on your computer, type:
#   python3 -m pip install name_of_library
 
import netCDF4           as nc     # To read NetCDF files.
import numpy             as np     # For mathematical operations on arrays.
import matplotlib.pyplot as plt    # To create figures.
from   matplotlib import colormaps # To use specific color palettes.
 
# Read from the atlas the different variables that we will need:
 
atla = nc.Dataset( 'atlas_chesbay_v20240319.nc' )
long = atla.variables['longitude'         ][:] # Shape: (         348)
lati = atla.variables['latitude'          ][:] # Shape: (    567     )
temp = atla.variables['temperature_bottom'][:] # Shape: (12, 567, 348) (deg.C)
omca = atla.variables['Omega_ca_bottom'   ][:] # Shape: (12, 567, 348)
del atla                                       # No longer necessary.
 
# Read the longitudes/latitudes corresponding to the COastLIne of the Bay:
 
clnc = nc.Dataset( 'coastlines_chesbay.nc' )
coli = clnc.variables['coli'][:] # Shape: (2, 314078) (longitudes, latitudes)
del clnc                         # No longer necessary.
 
# Parameters for equation of net calcification rate (Rivest & Brush, pers.comm.)
 
a    = 158.20
b    = 1.0520
c    = 0.9323
 
# NET Calcification rate (netc), units: g CaCO3 / (g tissue DW) / day.
 
netc = 1. - np.exp( - b * (omca - c) )
netc = netc * a * 1.e-3 * np.exp( 0.0271 * (temp - 25.) )
del a,b,c,omca,temp # No longer necessary.
 
# Produce two FIGUres, the first for January, second for September.
 
for figu in range(1, 3):
  if   figu == 1:
    mont = 0      # FIGUre 1 => MONTh of January   (index 0).
  elif figu == 2:
    mont = 8      # FIGUre 2 => MONTh of September (index 8).
 
# Define a figure and its axes:
 
  fig, axis = plt.subplots()
 
# Plot the calcification rate: specify the lon/lat of the atlas (long,lati),
#                              select the current MONTh,
#                              use a divergent ColorMAP (ReD to BlUe, Reversed),
#                              set the MIN/MAX Values of the figure to +/- 0.2.
 
  hand = axis.pcolormesh( long, lati, netc[mont,:,:], cmap = colormaps['RdBu_r'], vmin = - 0.2, vmax = 0.2 )
 
# Label the month on the top-right of the figure:
 
  if   mont == 0:
    axis.text( - 77.1, 39.4, 'January'   )
  elif mont == 8:
    axis.text( - 77.1, 39.4, 'September' )
 
# Add a colorbar with a suitable label:
 
  cbar = fig.colorbar( hand )
  cbar.set_label( 'Net calcification rate (g CaCO$_3$ (g tissue DW)$^{-1}$ day$^{-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:
 
  axis.set_xlim([ - 77.4, - 75.1 ]) # Longitudes (full extent of atlas)
  axis.set_ylim([   36.6,   39.6 ]) # Latitudes  (full extent of atlas)
# axis.set_xlim([ - 76.5, - 76.1 ]) # Longitudes (Elizabeth River     )
# axis.set_ylim([   36.7,   37.2 ]) # Latitudes  (Elizabeth River     )
 
# The 567 x 348 pixels of the atlas are approximately 600m x 600m in size.
# By enforcing a 'equal' ratio between the horizontal/vertical axes,
#   we are giving a reasonable aspect ratio to our map:
 
  plt.gca().set_aspect( 'equal', adjustable = 'box' )
 
# Overlay the COastLIne of the Chesapeake Bay on top of the calcification rate:

  plt.plot( coli[0,:], coli[1,:], color = 'gray', linewidth = 0.5 )

plt.show() # Show the two months on the screen for comparison.