{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"# Calculating Temporal Averages with GeoCAT-comp vs Xarray\n",
"\n",
"With temporally large datasets, computing seasonal and annual averages are a great ways to summarize the data and make it easier to manage and understand. You may want to take hourly, daily, or monthly data and compute seasonal or annual averages."
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Challenges\n",
"When using data that has a daily or finer resolution (e.g. hourly), calculating an annual average is simple. Every day and hour has the same length, so an unweighted average will work.\n",
"\n",
"But when using data that is monthly, things can get a bit tricky. Not every month is created equal. February has 28 or 29 days and March has 31 days. Since monthly data has one value for each month, those points can't be averaged in the usual way. A weighted average is needed.\n",
"\n",
"While it is tempting to quickly compute monthly to annual averages with `Xarray`'s `resample` or `groupby` functions, we need to be careful to specify the weights. Unfortunately, `Xarray` doesn't support weighted `resample` or `groupby` at the time this post was created, but `geocat-comp.climatology_average` builds upon `Xarray` to compute the weights for you.\n",
"\n",
"Below is a plot showing the difference between computing the winter average temperature from monthly data using the incorrect unweighted average and the correct weighted average.\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Demonstration\n",
"In this post, I'll show how to compute seasonal averages from monthly data the naive way (with unweighted averages) and the correct way (with weighted averages).\n",
"\n",
"### Imports"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"import cartopy as cart\n",
"import geocat.comp as gc\n",
"import geocat.viz as gv\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import xarray as xr"
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Helper function to make all of the plots the same way but with different data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"def custom_plot(data, title):\n",
" # Generate figure (set its size (width, height) in inches)\n",
" plt.figure(figsize=(14, 7))\n",
"\n",
" # Generate axes, using Cartopy\n",
" projection = cart.crs.PlateCarree()\n",
" ax = plt.axes(projection=projection)\n",
" ax.add_feature(cart.feature.LAND, zorder=10, edgecolor='k')\n",
"\n",
" # Draw coastlines\n",
" ax.coastlines()\n",
" ax.gridlines(alpha=0.5)\n",
"\n",
" if 'Difference' in title:\n",
" # Contourf-plot data (for filled contours)\n",
" p = data.plot.contourf(\n",
" ax=ax,\n",
" vmin=-0.1,\n",
" vmax=0.1,\n",
" levels=11,\n",
" cmap='bwr',\n",
" add_colorbar=False,\n",
" transform=projection,\n",
" extend='neither',\n",
" )\n",
"\n",
" # Add horizontal colorbar\n",
" cbar = plt.colorbar(p, orientation='horizontal', shrink=0.5)\n",
" cbar.ax.tick_params(labelsize=14)\n",
" cbar.set_ticks(np.linspace(-0.1, 0.1, 6))\n",
" else:\n",
" # Contourf-plot data (for filled contours)\n",
" p = data.plot.contourf(\n",
" ax=ax,\n",
" vmin=20,\n",
" vmax=30,\n",
" levels=11,\n",
" cmap='inferno',\n",
" add_colorbar=False,\n",
" transform=projection,\n",
" extend='neither',\n",
" )\n",
"\n",
" # Add horizontal colorbar\n",
" cbar = plt.colorbar(p, orientation='horizontal', shrink=0.5)\n",
" cbar.ax.tick_params(labelsize=14)\n",
" cbar.set_ticks(np.linspace(20, 30, 6))\n",
"\n",
" # Use geocat.viz.util convenience function to set axes tick values\n",
" gv.set_axes_limits_and_ticks(\n",
" ax,\n",
" xlim=(-180, -70),\n",
" ylim=(-20, 20),\n",
" xticks=np.arange(-180, -70, 10),\n",
" yticks=np.arange(-20, 20, 5),\n",
" )\n",
"\n",
" # Use geocat.viz.util convenience function to make plots look like NCL plots by using latitude, longitude tick labels\n",
" gv.add_lat_lon_ticklabels(ax)\n",
"\n",
" # Use geocat.viz.util convenience function to add minor and major tick lines\n",
" gv.add_major_minor_ticks(ax, labelsize=12)\n",
"\n",
" # Use geocat.viz.util convenience function to add titles to left and right of the plot axis.\n",
" gv.set_titles_and_labels(\n",
" ax,\n",
" maintitle=title,\n",
" lefttitle=\"Winter Average\",\n",
" lefttitlefontsize=16,\n",
" righttitle=data.units,\n",
" righttitlefontsize=16,\n",
" xlabel=\"\",\n",
" ylabel=\"\",\n",
" )\n",
"\n",
" # Show the plot\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Read in and format data\n",
"The data we will be using is a subset from RDA dataset ds277.0 - 'NOAA NCEP Optimum Interpolation Sea Surface Temperature Analysis'. It contains monthly average sea surface temperatures over the eastern equitorial Pacific from 1982 to 1986. We will be computing seasonal averages from this data and comparing the two different methods for doing this calculation."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"ds = xr.open_dataset('603321.sst.sst.mnmean.nc')\n",
"ds = ds.sst # Pull out the sea surface temperature data\n",
"ds = ds.isel(\n",
" time=range(1, 49)\n",
") # Remove the first data point so that we have an equal number of data points from each month"
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### The incorrect way to compute seasonal averages from monthly data"
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"It's easy to compute an unweighted average using `xarray` functionality; however, this generates inaccurate results. Here is what the ***incorrect*** way of doing this looks like. Notice that the result has a `season` dimension. These seasons are quarters of the year representing the meteorological seasons. (e.g. December, January, and February for winter)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/html": [
"