Climate variability

This second notebook showcases similar functionalities as example 1 but in addition it illustrates how RADWave can be used to estimate the influence of climate trends on regional wave patterns.

As in the previous notebook, we will be querying data already downloaded from the obtained from Australian Ocean Data Network portal AODN.

We recomend to look at RADWave documentation and the embeded video that explain how to select both a spatial bounding box and a temporal extent from the portal and how to export the file containing the List of URLs. This TXT file contains a list of NETCDF files for each available satellites.

Loading RADWave library and initialisation

We first start by importing RADwave library into our working space.

import io
import scipy
import requests
import numpy as np 
import pandas as pd

import RADWave as rwave

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

Here we will take the list of NetCDF URLs downloaded from the wave data portal containing the radar altimeter data for a region located offshore New Zealand [‘…/dataset/IMOS_NZ.txt’]

For a detail overview of the options available in the waveAnalysis class, you can have a look at the waveAnalysis API.

During initialisation we also specify:

  • bbox the bounding box of the geographical extent for the uploaded dataset following the convention [lon min,lon max,lat min,lat max]

  • stime the starting time of wave climate analysis following the convention [year, month, day] [we chose the 1st of January 1985]

  • etime the ending time of wave climate analysis following the convention [year, month, day] [we chose the 31st of December 2018]

Here again, we don’t specify a list of particular satellites to use (satNames keyword) so all of them will be queried. In other words we will look at all the records from the 10 altimeters:

JASON-2 - JASON-3 - SARAL - SENTINEL-3A - CRYOSAT-2 - ENVISAT - GEOSAT - ERS-2 - GFO - TOPEX.

wa = rwave.waveAnalysis(altimeterURL='../pracenv/dataset/IMOS_NZ.txt', bbox=[175.0,177.0,-47.0,-45.0], 
                  stime=[1985,1,1], etime=[2019,12,31])

Processing altimeters data

After class initialisation querying the actual dataset is realised by calling the processAltimeterData function. The description of this function is available from the API.

The function can take some times to execute depending on the number of NETCDF files to load and the size of the dataset to query (here it should not take more than 30 s).

RADWave uses the uploaded file containing the list of URLs to query via THREDDS the remote data. This operation can take several minutes and when looking at a large region it is recommended to divide the analyse in smaller regions and download a series of URLs text file instead of the entire domain directly.

wa.processAltimeterData(max_qc=1, altimeter_pick='all', saveCSV = 'altimeterData.csv')
Processing Altimeter Dataset 

   +  name JASON-2     / number of tracks                               4   
   +  name JASON-3     / number of tracks                               4   
   +  name SARAL       / number of tracks                               4   
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-3-a77cf9e50a80> in <module>
----> 1 wa.processAltimeterData(max_qc=1, altimeter_pick='all', saveCSV = 'altimeterData.csv')

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/RADWave/altiwave.py in processAltimeterData(self, max_qc, altimeter_pick, saveCSV)
    366                             time_m = dataframe.groupby(["year", "month", "day"])[
    367                                 ["time"]
--> 368                             ].apply(np.median)
    369                             time_m.name = "time"
    370                             qc_m = dataframe.groupby(["year", "month", "day"])[

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/pandas/core/groupby/groupby.py in apply(self, func, *args, **kwargs)
    857         with option_context("mode.chained_assignment", None):
    858             try:
--> 859                 result = self._python_apply_general(f, self._selected_obj)
    860             except TypeError:
    861                 # gh-20949

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/pandas/core/groupby/groupby.py in _python_apply_general(self, f, data)
    890             data after applying f
    891         """
--> 892         keys, values, mutated = self.grouper.apply(f, data, self.axis)
    893 
    894         return self._wrap_applied_output(

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/pandas/core/groupby/ops.py in apply(self, f, data, axis)
    184         ):
    185             try:
--> 186                 result_values, mutated = splitter.fast_apply(f, sdata, group_keys)
    187 
    188             except libreduction.InvalidApply as err:

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/pandas/core/groupby/ops.py in fast_apply(self, f, sdata, names)
    971         # must return keys::list, values::list, mutated::bool
    972         starts, ends = lib.generate_slices(self.slabels, self.ngroups)
--> 973         return libreduction.apply_frame_axis0(sdata, f, names, starts, ends)
    974 
    975     def _chop(self, sdata: DataFrame, slice_obj: slice) -> DataFrame:

pandas/_libs/reduction.pyx in pandas._libs.reduction.apply_frame_axis0()

<__array_function__ internals> in median(*args, **kwargs)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/lib/function_base.py in median(a, axis, out, overwrite_input, keepdims)
   3519     """
   3520     r, k = _ureduce(a, func=_median, axis=axis, out=out,
-> 3521                     overwrite_input=overwrite_input)
   3522     if keepdims:
   3523         return r.reshape(k)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/lib/function_base.py in _ureduce(a, func, **kwargs)
   3427         keepdim = (1,) * a.ndim
   3428 
-> 3429     r = func(a, **kwargs)
   3430     return r, keepdim
   3431 

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/lib/function_base.py in _median(a, axis, out, overwrite_input)
   3553             part = a
   3554     else:
-> 3555         part = partition(a, kth, axis=axis)
   3556 
   3557     if part.shape == ():

<__array_function__ internals> in partition(*args, **kwargs)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/core/fromnumeric.py in partition(a, kth, axis, kind, order)
    746     else:
    747         a = asanyarray(a).copy(order="K")
--> 748     a.partition(kth, axis=axis, kind=kind, order=order)
    749     return a
    750 

KeyboardInterrupt: 

Once the dataset has been queried, we can plot the altimeter data points on a map using the visualiseData function.

This function plots and saves in a figure the geographical coordinates of processed altimeter data.

wa.visualiseData(title="Altimeter data tracks offshore New Zealand", extent=[166.,179.,-48.,-42.], 
                 addcity=None, markersize=30, zoom=8,
                 fsize=(10, 9), fsave=None)
/Users/getafix/anaconda3/envs/coast/lib/python3.6/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/Users/getafix/anaconda3/envs/coast/lib/python3.6/site-packages/cartopy/mpl/gridliner.py:331: UserWarning: The .ylabels_left attribute is deprecated. Please use .left_labels to toggle visibility instead.
  warnings.warn('The .ylabels_left attribute is deprecated. Please '
../_images/climatevariability_7_1.svg

Computing wave regime and seasonality

Time series

To perform wave analysis and compute the wave parameters discussed in the documentation, we run the generateTimeSeries function. And we then plot these time series by calling the plotTimeSeries function.

timeseries = wa.generateTimeSeries()
wa.plotTimeSeries(time=[1995,2019], series='H', fsize=(12, 5), fsave=None)
wa.plotTimeSeries(time=[1995,2019], series='T', fsize=(12, 5), fsave=None)
Max wave height: 8.944 m
Mean wave height: 3.212 m
Median wave height: 2.985 m
95th percentile wave height: 5.532 m
../_images/climatevariability_10_1.svg
Max wave period: 12.887 s
Mean wave period: 7.382 s
Median wave period: 7.275 s
95th percentile wave period: 9.440 s
../_images/climatevariability_10_3.svg

20-years analysis of the impact of climate trend

Oscillation in atmospheric patterns is known to alter regional weather conditions and associated trends in wave climate [Godoi et al., 2016].

Here we illustrate how the results obtained with RADWave can be used to investigate how climate patterns may affect wave parameters.

For the sake of the demonstration, we will focus our analysis on the following indices:

We first load the data associated to each index using Pandas functionalities.

# Defining the 20-years timeframe
time = [1998,2018] 

Monthly means of the SOI, AOI & SAMI index are sourced from the National Oceanic and Atmospheric Administration (NOAA) and the Natural Environment Research Council from the British Antarctic Survey (NERC).

For each climate index the anomalies are computed by subtracting overall mean from the monthly means.

Then, the same is done for the wave parameters in order to investigate how they are modulated by the climate modes.

SOI - Southern Oscillation Index

# Dataset URL
url = "https://www.ncdc.noaa.gov/teleconnections/enso/indicators/soi/data.csv"

# Using Pandas to load the file content
soi = requests.get(url).content
soi_data = pd.read_csv(io.StringIO(soi.decode('utf-8')),skiprows=1)

# Define year and month of each record
soi_data['year'] = soi_data['Date'] // 100
soi_data['month'] = soi_data['Date'] % 100 

# Extract the information for the specified time interval 
soi_df = soi_data.drop(soi_data[soi_data.year < time[0]].index)
soi_df = soi_df.drop(soi_df[soi_df.year > time[1]].index)

# Calculate the 20-years mean 
soi_mean = soi_df['Value'].mean()

# Compute and store the anomalies in the dataframe
soi_df['anomaly'] = soi_df['Value']-soi_mean

AOI - Antarctic Oscillation Index

# Dataset URL
url = "https://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/aao/monthly.aao.index.b79.current.ascii"

# Using Pandas to load the file content
aoi = requests.get(url).content
aoi_data = pd.read_csv(io.StringIO(aoi.decode('utf-8')),delimiter=r"\s+",header=None)

# Rename columns to fit with RADWave dataframe
aoi_data = aoi_data.rename(columns={0:"year", 1:"month", 2:"Value"})

# Extract the information for the specified time interval 
aoi_df = aoi_data.drop(aoi_data[aoi_data.year < time[0]].index)
aoi_df = aoi_df.drop(aoi_df[aoi_df.year > time[1]].index)

# Calculate the 20-years mean 
aoi_mean = aoi_df['Value'].mean()

# Compute and store the anomalies in the dataframe
aoi_df['anomaly'] = aoi_df['Value']-aoi_mean

SAMI - Southern Annular Mode Index

# Dataset URL
url = "http://www.nerc-bas.ac.uk/public/icd/gjma/newsam.1957.2007.txt"

# Using Pandas to load the file content
sam = requests.get(url).content
sam_data = pd.read_csv(io.StringIO(sam.decode('utf-8')),delimiter=r"\s+")

# Rename month values to fit with RADWave dataframe
sam_data = sam_data.rename(columns={"JAN":1, "FEB":2,
                      "MAR":3, "APR":4,
                      "MAY":5, "JUN":6,
                      "JUL":7, "AUG":8,
                      "SEP": 9, "OCT":10,
                      "NOV":11, "DEC":12})

# Rename columns to fit with RADWave dataframe
sam_data = sam_data.stack().reset_index()
sam_data = sam_data.rename(columns={"level_0":"year", "level_1":"month", 0:"Value"})


# Extract the information for the specified time interval 
sam_df = sam_data.drop(sam_data[sam_data.year < time[0]].index)
sam_df = sam_df.drop(sam_df[sam_df.year > time[1]].index)

# Calculate the 20-years mean 
sam_mean = sam_df['Value'].mean()

# Compute and store the anomalies in the dataframe
sam_df['anomaly'] = sam_df['Value']-sam_mean

RADWave significant wave height and wave period anomalies

Significant wave height

# Get monthly significant wave height 
wh_data = wa.timeseries.groupby(['year', 'month'])[['wh']].apply(np.mean).reset_index()

# Extract the information for the specified time interval 
wh_df = wh_data.drop(wh_data[wh_data.year < time[0]].index)
wh_df = wh_df.drop(wh_df[wh_df.year > time[1]].index)

# Calculate the 20-years mean 
wh_mean = wh_df['wh'].mean()

# Compute and store the anomalies in the dataframe
wh_df['anomaly'] = wh_df['wh']-wh_mean

Wave period

# Get monthly mean wave period 
T_data = wa.timeseries.groupby(['year', 'month'])[['period']].apply(np.mean).reset_index()

# Extract the information for the specified time interval 
T_df = T_data.drop(T_data[T_data.year < time[0]].index)
T_df = T_df.drop(T_df[T_df.year > time[1]].index)

# Calculate the 20-years mean 
T_mean = T_df['period'].mean()

# Compute and store the anomalies in the dataframe
T_df['anomaly'] = T_df['period']-T_mean

Correlations

Monthly mean anomalies of significant wave height and wave period can be correlated with monthly mean anomaly time series of the SOI, AOI & SAMI index by computing the Pearson’s correlation coefficient ® for the region of interest.

We use scipy.stats.pearsonr function to make this calculation. This function returns 2 values:

  • r: Pearson’s correlation coefficient.

  • p: Two-tailed p-value.

Examples of Pearson’s correlation coefficient calculation between climate indices and the significant wave height are provided below:

# Pearson correlation between significant wave height and SOI
monthly_wh_soi = scipy.stats.pearsonr(soi_df['anomaly'],wh_df['anomaly']) 
print('+ Pearson correlation between significant wave height and SOI:',monthly_wh_soi[0],'\n')

# Pearson correlation between significant wave height and AOI
monthly_wh_aoi = scipy.stats.pearsonr(aoi_df['anomaly'],wh_df['anomaly']) 
print('+ Pearson correlation between significant wave height and AOI:',monthly_wh_aoi[0],'\n')

# Pearson correlation between significant wave height and SAMI
monthly_wh_sam = scipy.stats.pearsonr(sam_df['anomaly'],wh_df['anomaly']) 
print('+ Pearson correlation between significant wave height and SAMI:',monthly_wh_sam[0],'\n')
+ Pearson correlation between significant wave height and SOI: -0.08908662862356329 

+ Pearson correlation between significant wave height and AOI: -0.13649468119337463 

+ Pearson correlation between significant wave height and SAMI: -0.011441259490892312