eReefs dataset cross section¶
We have seen how to make a map of the eReefs dataset. In some occasions, we might be interested in looking at the evolution of a specific parameter over a cross-section. Here we will see how this can be simply done with the librairies we have seen so far…
Load the required Python libraries¶
First of all, load the necessary libraries:
numpy
matplotlib
cartopy
import os
import numpy as np
import datetime as dt
import netCDF4
from netCDF4 import Dataset, num2date
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
cartopy.config['data_dir'] = os.getenv('CARTOPY_DIR', cartopy.config.get('data_dir'))
import cmocean
from matplotlib import pyplot as plt
# %config InlineBackend.figure_format = 'retina'
plt.ion() # To trigger the interactive inline mode
Connect to the OPeNDAP endpoint for a specified month.¶
We query the server based on the OPeNDAP URL for a specific file. We will use the data from the AIMS server.
gbr4: we use the 4km model and daily data for a month specified
month = 4
year = 2019
netCDF_datestr = str(year)+'-'+format(month, '02')
print('File chosen time interval:',netCDF_datestr)
# GBR4 HYDRO
#inputFile = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-"+netCDF_datestr+".nc"
# GBR4 BIO
inputFile = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/GBR4_H2p0_B3p1_Cq3b_Dhnd/daily-monthly/EREEFS_AIMS-CSIRO_GBR4_H2p0_B3p1_Cq3b_Dhnd_bgc_daily-monthly-"+netCDF_datestr+".nc"
File chosen time interval: 2019-04
We now load the dataset within the Jupyter environment…
nc_data = Dataset(inputFile, 'r')
print('Get the list of variable in the file:')
print(list(nc_data.variables.keys()))
ncdata = nc_data.variables
Get the list of variable in the file:
['alk', 'BOD', 'Chl_a_sum', 'CO32', 'DIC', 'DIN', 'DIP', 'DOR_C', 'DOR_N', 'DOR_P', 'Dust', 'EFI', 'FineSed', 'Fluorescence', 'HCO3', 'Kd_490', 'MPB_Chl', 'MPB_N', 'Mud-carbonate', 'Mud-mineral', 'Nfix', 'NH4', 'NO3', 'omega_ar', 'Oxy_sat', 'Oxygen', 'P_Prod', 'PAR', 'PAR_z', 'pco2surf', 'PH', 'PhyL_Chl', 'PhyL_N', 'PhyS_Chl', 'PhyS_N', 'PhyS_NR', 'PIP', 'salt', 'TC', 'temp', 'TN', 'TP', 'Tricho_Chl', 'Tricho_N', 'Z_grazing', 'ZooL_N', 'ZooS_N', 'zc', 'time', 'latitude', 'longitude', 'CH_N', 'CS_bleach', 'CS_Chl', 'CS_N', 'EpiPAR_sg', 'eta', 'MA_N', 'MA_N_pr', 'month_EpiPAR_sg', 'R_400', 'R_410', 'R_412', 'R_443', 'R_470', 'R_486', 'R_488', 'R_490', 'R_510', 'R_531', 'R_547', 'R_551', 'R_555', 'R_560', 'R_590', 'R_620', 'R_640', 'R_645', 'R_665', 'R_667', 'R_671', 'R_673', 'R_678', 'R_681', 'R_709', 'R_745', 'R_748', 'R_754', 'R_761', 'R_764', 'R_767', 'R_778', 'Secchi', 'SG_N', 'SG_N_pr', 'SG_shear_mort', 'SGD_N', 'SGD_N_pr', 'SGD_shear_mort', 'SGH_N', 'SGH_N_pr', 'SGHROOT_N', 'SGROOT_N', 'TSSM', 'Zenith2D']
To get information for a specific variables, we can use the following:
ncdata['alk']
<class 'netCDF4._netCDF4.Variable'>
float32 alk(time, k, latitude, longitude)
coordinates: time zc latitude longitude
short_name: alk
units: mmol m-3
long_name: Total alkalinity
_ChunkSizes: [ 1 1 133 491]
unlimited dimensions: time
current shape = (30, 17, 723, 491)
filling off
ncdata['alk'].long_name
'Total alkalinity'
ncdata['alk'].units
'mmol m-3'
Load variables¶
We then load the longitude/latidude in the Jupyter environment:
# Starting with the spatial domain
lat = ncdata['latitude'][:]
lon = ncdata['longitude'][:]
print('eReefs model spatial extent:\n')
print(' - Longitudinal extent:',lon.min(),lon.max())
print(' - Latitudinal extent:',lat.min(),lat.max())
eReefs model spatial extent:
- Longitudinal extent: 142.168788 156.868788
- Latitudinal extent: -28.696022 -7.036022
We can also check the temporal and z-coordinate extent of the chosem dataset:
# Get time span of the dataset
time_var = ncdata['time']
# Starting time
dtime = netCDF4.num2date(time_var[0],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print(' - start time: ',daystr,'\n')
# Ending time
dtime = netCDF4.num2date(time_var[-1],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print(' - end time: ',daystr,'\n')
ntime = len(time_var)
print(' - Number of time steps',ntime,'\n')
zc = ncdata['zc'][:]
nlay = len(zc)
print(' - Number of vertical layers',nlay)
- start time: 2019-Apr-01 02:00
- end time: 2019-Apr-30 02:00
- Number of time steps 30
- Number of vertical layers 17
Cross-section¶
Let’s make 2 cross-sections along specific latitude and longitude coordinates.
chosen latitude: 11 S
chosen longitude: 148 E
The AIMS dataset has been resample on a regular grid, so first we find the closest point to the desired latitude.
def find_nearest(array, value):
'''
Find index of nearest value in a numpy array
'''
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
We will then extract all the point along this latitude:
def eReefs_cross_data(data, tstep, latID=None, lonID=None):
'''
Extract specified dataset along a given latitude or longitude.
args:
- dataname: specified variable name
- latID: latitudinal index
- lonID: longitudinal index
'''
# Get data
if latID is not None:
return data[tstep, :, latID, :]
elif lonID is not None:
return data[tstep, :, :, lonID]
Now we will call these functions…
selectedLatIndex = find_nearest(lat, -11.)
selectedLonIndex = find_nearest(lon, 148)
selectedTimeIndex = 29
selectedVariable = 'alk'
alkLat = eReefs_cross_data(ncdata[selectedVariable], selectedTimeIndex,
selectedLatIndex, None)
alkLon = eReefs_cross_data(ncdata[selectedVariable], selectedTimeIndex,
None, selectedLonIndex)
print(' ')
print('Data range: ')
print(np.nanmin(alkLat),np.nanmax(alkLat))
Data range:
2183.413 2270.2312
We then define the plotting function:
def plot_cross_section(ncdata, tstep, dataname, xdata, ydata, color, size, fname,
xext=None, dext=None, show=False, save=False):
'''
This function plots a cross-section for a specified time interval along a given latitude or longitude.
args:
- ncdata: netcdf dataset
- tstep: specified time index
- dataname: specified variable name
- xdata: x-axis along the desired cross-section
- ydata: data along the desired cross-section
- xext: extent to plot for the x-axis
- dext: data minimum and maximum values for plotting
- color: colormap to use for the plot (here one can use the cmocean library: https://matplotlib.org/cmocean/#installation)
- size: figure size
- fname: figure name when saved on disk, it is worth noting that the specified time index will be automatically added
- show: set to True when the map is shown in the jupyter environment directly
- save: set to True to save figure on disk
'''
fig = plt.figure(figsize=size, facecolor='w', edgecolor='k')
ax = plt.axes()
plt.xlabel('Cross-section (degree)')
plt.ylabel('Depth (m)')
if xext is not None:
plt.xlim(xext[0], xext[1])
zc = ncdata['zc'][:]
plt.ylim(zc.min(), zc.max())
if dext is not None:
cm = plt.pcolormesh(xdata, zc, ydata[:,:],
cmap = color,
vmin = dext[0],
vmax = dext[1],
edgecolors = 'face',
shading ='auto')
else:
cm = plt.pcolormesh(xdata, zc, ydata[:,:],
cmap = color,
edgecolors = 'face',
shading ='auto')
dtime = netCDF4.num2date(ncdata['time'][tstep],ncdata['time'].units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
plt.title(ncdata[dataname].long_name+', %s UTC+10' % (daystr), fontsize=11);
# Color bar
cbar = fig.colorbar(cm, ax=ax, fraction=0.01, pad=0.025)
cbar.set_label(ncdata[dataname].units, rotation=90, labelpad=5, fontsize=10)
cbar.ax.tick_params(labelsize=8)
# Get z-coordinate lines
for k in range(len(zc)):
plt.plot(xext,[zc[k],zc[k]],lw=0.5,c='k',alpha=0.25)
if show:
if save:
plt.savefig(f"{fname}_cross_time{tstep:04}.png",dpi=300,
bbox_inches='tight')
plt.tight_layout()
plt.show()
else:
plt.savefig(f"{fname}_cross_time{tstep:04}.png",dpi=300,
bbox_inches='tight')
fig.clear()
plt.close(fig)
plt.clf()
return
color = cmocean.cm.matter
fname = 'lat11'
xext = [lon.min(),150.25]
size = (9,3)
dext = [2190.,2270.]
plot_cross_section(ncdata, selectedTimeIndex, selectedVariable, lon, alkLat, color, size,
fname, xext, dext, show=True, save=False)
fname = 'lon148'
xext = [-20.5,-11.5]
size = (10,3)
plot_cross_section(ncdata, selectedTimeIndex, selectedVariable, lat, alkLon, color, size,
fname, xext, dext, show=True, save=False)
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
Build a function to simplify the process¶
def eReefs_cross(ncdata, tstep, dataname, color, latVal=None, lonVal=None, size=None,
fname=None, xext=None, dext=None, show=False, save=False):
'''
This function plots cross-sections for a specified time interval along a given latitude and/or longitude.
args:
- ncdata: netcdf dataset
- tstep: specified time index
- dataname: specified variable name
- color: colormap to use for the plot (here one can use the cmocean library: https://matplotlib.org/cmocean/#installation)
- latVal: latitude of the desired cross-section
- lonVal: longitude of the desired cross-section
- size: figure size
- fname: figure name when saved on disk, the time step value is automatically added
- xext: extent to plot for the x-axis for the cross-section
- dext: data minimum and maximum values for plotting
- show: set to True when the map is shown in the jupyter environment directly
- save: set to True to save figure on disk
'''
lat = ncdata['latitude'][:]
lon = ncdata['longitude'][:]
if latVal is not None:
LatID = find_nearest(lat, latVal)
dataLat = eReefs_cross_data(ncdata[dataname], tstep, LatID, None)
plot_cross_section(ncdata, tstep, dataname, lon, dataLat, color, size,
fname, xext, dext, show, save)
if lonVal is not None:
LonID = find_nearest(lon, lonVal)
dataLon = eReefs_cross_data(ncdata[dataname], tstep, None, LonID)
plot_cross_section(ncdata, tstep, dataname, lat, dataLon, color, size,
fname, xext, dext, show, save)
return
selectedTimeIndex = 10
selectedVariable = 'alk'
color = cmocean.cm.matter
latVal = -11.
lonVal = 148.
xLon = [142,150.25]
xLat = [-20.5,-11.5]
dext = [2190.,2270.]
fLat = 'lat11'
fLon = 'lon148'
sizeLat = (9,3)
sizeLon = (10,3)
eReefs_cross(ncdata, selectedTimeIndex, selectedVariable, color,
latVal, None, sizeLat, fLat, xLon, dext, show=True,
save=False)
eReefs_cross(ncdata, selectedTimeIndex, selectedVariable, color,
None, lonVal, sizeLon, fLon, xLat, dext, show=True,
save=False)
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>