03-summarize-data-seasons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"earth-lab-logo-rgb.png\" width=\"150\" height=\"150\" />\n",
    "\n",
    "# Earth Analytics Education - Climate 101 Workshop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <i class=\"fa fa-ship\" aria-hidden=\"true\"></i> Chapter 13 - NETCDF 4 Climate Data in Open Source Python \n",
    "\n",
    "In this chapter, you will learn how to work with Climate Data Sets (MACA v2 for the United states) stored in netcdf 4 format using open source **Python**.\n",
    "\n",
    "## <i class=\"fa fa-graduation-cap\" aria-hidden=\"true\"></i> Learning Objectives\n",
    "\n",
    "After completing this chapter, you will be able to:\n",
    "\n",
    "* Summarize MACA v 2 climate data stored in netcdf 4 format by seasons across all time periods using `xarray`.\n",
    "* Summarize MACA v 2 climate data stored in netcdf 4 format by seasons and across years using `xarray`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Seasonal Averages Using MACA vs Climate Data \n",
    "\n",
    "In this lesson, you will learn how to calculate seasonal averages over several years\n",
    "using MACA v 2 Climate Data downloaded in `netcdf4` format using \n",
    "`xarray`.\n",
    "\n",
    "In this example you will use the forecast temperature data downloaded from the \n",
    "northwestknowledge.net website. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure ea dir exists - if not create it\n",
    "import os\n",
    "import earthpy as et\n",
    "ea_path = os.path.join(et.io.HOME, 'earth-analytics', 'data')\n",
    "\n",
    "if not os.path.exists(ea_path):\n",
    "    os.makedirs(ea_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import seaborn as sns\n",
    "import geopandas as gpd\n",
    "import earthpy as et\n",
    "import xarray as xr\n",
    "import regionmask\n",
    "\n",
    "# Plotting options\n",
    "sns.set(font_scale=1.3)\n",
    "sns.set_style(\"white\")\n",
    "\n",
    "# Optional - set your working directory if you wish to use the data\n",
    "# accessed lower down in this notebook (the USA state boundary data)\n",
    "os.chdir(os.path.join(et.io.HOME,\n",
    "                      'earth-analytics',\n",
    "                      'data'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin, you can download and open up a MACA v2 netcdf file. The file below is a \n",
    "projected maximum temperature dataset downscaled using the `BNU-ESM` model for 2006-2099."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get netcdf file\n",
    "data_path_monthly = 'http://thredds.northwestknowledge.net:8080/thredds/dodsC/agg_macav2metdata_tasmax_BNU-ESM_r1i1p1_rcp45_2006_2099_CONUS_monthly.nc'\n",
    "\n",
    "# Open up the data\n",
    "with xr.open_dataset(data_path_monthly) as file_nc:\n",
    "    monthly_forecast_temp_xr = file_nc\n",
    "\n",
    "# xarray object\n",
    "monthly_forecast_temp_xr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example below you subset data for the state of \n",
    "California similar to what you did in the previous lesson. You can \n",
    "select any state that you wish for this analysis!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download natural earth data to generate AOI\n",
    "et.data.get_data(\n",
    "    url=\"https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_1_states_provinces_lakes.zip\")\n",
    "\n",
    "states_path = \"earthpy-downloads/ne_50m_admin_1_states_provinces_lakes\"\n",
    "states_path = os.path.join(\n",
    "    states_path, \"ne_50m_admin_1_states_provinces_lakes.shp\")\n",
    "\n",
    "states_gdf = gpd.read_file(states_path)\n",
    "cali_aoi = states_gdf[states_gdf.name == \"California\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Would this be better if it only returned 4 values?? probably so\n",
    "# Helper Function to extract AOI\n",
    "def get_aoi(shp, world=True):\n",
    "    \"\"\"Takes a geopandas object and converts it to a lat/ lon\n",
    "    extent \"\"\"\n",
    "\n",
    "    lon_lat = {}\n",
    "    # Get lat min, max\n",
    "    aoi_lat = [float(shp.total_bounds[1]), float(shp.total_bounds[3])]\n",
    "    aoi_lon = [float(shp.total_bounds[0]), float(shp.total_bounds[2])]\n",
    "\n",
    "    # Handle the 0-360 lon values\n",
    "    if world:\n",
    "        aoi_lon[0] = aoi_lon[0] + 360\n",
    "        aoi_lon[1] = aoi_lon[1] + 360\n",
    "    lon_lat[\"lon\"] = aoi_lon\n",
    "    lon_lat[\"lat\"] = aoi_lat\n",
    "    return lon_lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get lat min, max from Cali aoi extent\n",
    "cali_bounds = get_aoi(cali_aoi)\n",
    "\n",
    "# Slice by time & aoi location\n",
    "start_date = \"2059-12-15\"\n",
    "end_date = \"2099-12-15\"\n",
    "\n",
    "cali_temp = monthly_forecast_temp_xr[\"air_temperature\"].sel(\n",
    "    time=slice(start_date, end_date),\n",
    "    lon=slice(cali_bounds[\"lon\"][0], cali_bounds[\"lon\"][1]),\n",
    "    lat=slice(cali_bounds[\"lat\"][0], cali_bounds[\"lat\"][1]))\n",
    "cali_temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Time Period start: \", cali_temp.time.min().values)\n",
    "print(\"Time Period end: \", cali_temp.time.max().values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the region mask object - this is used to identify each region\n",
    "# cali_region = regionmask.from_geopandas(cali_aoi,\n",
    "#                                         names=\"name\",\n",
    "#                                         name=\"name\",\n",
    "#                                         abbrevs=\"iso_3166_2\")\n",
    "cali_mask = regionmask.mask_3D_geopandas(cali_aoi,\n",
    "                                         monthly_forecast_temp_xr.lon,\n",
    "                                         monthly_forecast_temp_xr.lat)\n",
    "# Mask the netcdf data\n",
    "cali_temp_masked = cali_temp.where(cali_mask)\n",
    "cali_temp_masked.dims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cali_temp.values.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the mean temperature for each season across the entire dataset. \n",
    "This will produce 4 arrays - one representing mean temperature for each seasons. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cali_season_summary = cali_temp_masked.groupby(\n",
    "    'time.season').mean('time', skipna=True)\n",
    "\n",
    "# This will create 4 arrays - one for each season showing mean temperature values\n",
    "cali_season_summary.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the seasonal data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "caption": "Facet plot showing seasonal mean temperature values for the State of California summarized over time."
   },
   "outputs": [],
   "source": [
    "# Create a plot showing mean temperature aross seasons\n",
    "cali_season_summary.plot(col='season', col_wrap=2, figsize=(10, 10))\n",
    "plt.suptitle(\"Mean Temperature Across All Selected Years By Season \\n California, USA\",\n",
    "             y=1.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate UnWeighted Seasonal Averages For By Season Across Each Year\n",
    "\n",
    "Above you created one single value per season which summarized seasonal data across all years.\n",
    "However you may want to look at seasonal variation year to year in the projected data.\n",
    "You can calculate seasonal statistcs by \n",
    "\n",
    "1. resampling the data and then\n",
    "2. grouping the data and summarizing it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resample the data by season across all years\n",
    "cali_season_mean_all_years = cali_temp_masked.resample(\n",
    "    time='QS-DEC', keep_attrs=True).mean()\n",
    "cali_season_mean_all_years.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summarize each array into one single (mean) value\n",
    "cali_seasonal_mean = cali_season_mean_all_years.groupby('time').mean([\n",
    "    \"lat\", \"lon\"])\n",
    "cali_seasonal_mean.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This data now has one value per season rather than an array\n",
    "cali_seasonal_mean.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "caption": "Scatter plot showing seasonal mean temperature values for the State of California displayed over time."
   },
   "outputs": [],
   "source": [
    "# Plot the data\n",
    "f, ax = plt.subplots(figsize=(10, 4))\n",
    "cali_seasonal_mean.plot(marker=\"o\",\n",
    "                        color=\"grey\",\n",
    "                        markerfacecolor=\"purple\",\n",
    "                        markeredgecolor=\"purple\")\n",
    "ax.set(title=\"Seasonal Mean Temperature\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export Seasonal Climate Project Data To .csv File\n",
    "At this point you can convert the data to a dataframe and export \n",
    "it to a `.csv` format if you wish."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to a dataframe\n",
    "cali_seasonal_mean_df = cali_seasonal_mean.to_dataframe()\n",
    "cali_seasonal_mean_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export a csv file\n",
    "cali_seasonal_mean_df.to_csv(\"cali-seasonal-temp.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Seasonal Data By Season \n",
    "\n",
    "Using `groupby()` you can group the data and plot it by season to better look at seasonal trends. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "caption": "Scatter plot showing seasonal mean temperature values for the State of California summarized over time and colored by season."
   },
   "outputs": [],
   "source": [
    "colors = {3: \"grey\", 6: \"lightgreen\", 9: \"green\", 12: \"purple\"}\n",
    "seasons = {3: \"spring\", 6: \"summer\", 9: \"fall\", 12: \"winter\"}\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10, 7))\n",
    "for month, arr in cali_seasonal_mean.groupby('time.month'):\n",
    "    arr.plot(ax=ax,\n",
    "             color=\"grey\",\n",
    "             marker=\"o\",\n",
    "             markerfacecolor=colors[month],\n",
    "             markeredgecolor=colors[month],\n",
    "             label=seasons[month])\n",
    "\n",
    "ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "ax.set(title=\"Seasonal Change in Mean Temperature Over Time\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weighted Summary by Season\n",
    "\n",
    "To begin, you will generate a list of days in each month which will be used to weight\n",
    "your seasonal summary data according to the the days in each month. \n",
    "\n",
    "* TODO -- redo this section to use the approach above which is perfect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate seasonal averages\n",
    "# http://xarray.pydata.org/en/stable/examples/monthly-means.html\n",
    "\n",
    "month_length = cali_temp_masked.time.dt.days_in_month\n",
    "month_length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, divide the data grouped by season by the total days represented\n",
    "in each season to create weighted values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "# https://stackoverflow.com/questions/59234745/is-there-any-easy-where-way-to-compute-seasonal-mean-with-xarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is returning values of 0 rather than na\n",
    "# Calculate a weighted mean by season\n",
    "cali_weighted_mean = ((cali_temp * month_length).resample(time='QS-DEC').sum() /\n",
    "                      month_length.resample(time='QS-DEC').sum())\n",
    "# Replace 0 values with nan\n",
    "cali_weighted_mean = cali_weighted_mean.where(cali_weighted_mean)\n",
    "cali_weighted_mean.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cali_weighted_season_value = cali_weighted_mean.groupby('time').mean([\n",
    "    \"lat\", \"lon\"])\n",
    "cali_weighted_season_value.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "caption": "Comparison of seasonal mean values - weighted vs unweighted by days in each month and colored by season."
   },
   "outputs": [],
   "source": [
    "colors = {3: \"grey\", 6: \"lightgreen\", 9: \"green\", 12: \"purple\"}\n",
    "seasons = {3: \"Spring\", 6: \"Summer\", 9: \"Fall\", 12: \"Winter\"}\n",
    "\n",
    "f, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharey=True)\n",
    "for month, arr in cali_weighted_season_value.groupby('time.month'):\n",
    "    arr.plot(ax=ax1,\n",
    "             color=\"grey\",\n",
    "             marker=\"o\",\n",
    "             markerfacecolor=colors[month],\n",
    "             markeredgecolor=colors[month],\n",
    "             label=seasons[month])\n",
    "\n",
    "ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "ax1.set(title=\"Weighted Seasonal Change in Mean Temperature Over Time\")\n",
    "\n",
    "for month, arr in cali_seasonal_mean.groupby('time.month'):\n",
    "    arr.plot(ax=ax2,\n",
    "             color=\"grey\",\n",
    "             marker=\"o\",\n",
    "             markerfacecolor=colors[month],\n",
    "             markeredgecolor=colors[month],\n",
    "             label=seasons[month])\n",
    "\n",
    "ax2.set(title=\"Unweighted Seasonal Change in Mean Temperature Over Time\")\n",
    "f.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want, you can compare the difference between weighted vs unweighted \n",
    "values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What does the difference look like weighted vs unweighted?\n",
    "cali_seasonal_mean - cali_weighted_season_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Same Analysis for the West Coast\n",
    "\n",
    "Above you calculate seasonal summaries for the state of California. You can \n",
    "implement the same analysis for each aoi region in a shapefile if you want \n",
    "following the workflow that you learned in the previous lesson. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create AOI Subset\n",
    "cali_or_wash_nev = states_gdf[states_gdf.name.isin(\n",
    "    [\"California\", \"Oregon\", \"Washington\", \"Nevada\"])]\n",
    "west_bounds = get_aoi(cali_or_wash_nev)\n",
    "\n",
    "# Create the mask\n",
    "west_mask = regionmask.mask_3D_geopandas(cali_or_wash_nev,\n",
    "                                         monthly_forecast_temp_xr.lon,\n",
    "                                         monthly_forecast_temp_xr.lat)\n",
    "\n",
    "\n",
    "# Slice by time & aoi location\n",
    "start_date = \"2059-12-15\"\n",
    "end_date = \"2099-12-15\"\n",
    "\n",
    "west_temp = monthly_forecast_temp_xr[\"air_temperature\"].sel(\n",
    "    time=slice(start_date, end_date),\n",
    "    lon=slice(west_bounds[\"lon\"][0], west_bounds[\"lon\"][1]),\n",
    "    lat=slice(west_bounds[\"lat\"][0], west_bounds[\"lat\"][1]))\n",
    "\n",
    "# Apply the mask\n",
    "west_temp_masked = west_temp.where(west_mask)\n",
    "west_temp_masked\n",
    "# Resample the data by season across all years\n",
    "#west_season_mean_all_years = west_temp_masked.groupby('region').resample(time='QS-DEC', keep_attrs=True).mean()\n",
    "# cali_seasonal_mean = cali_season_mean_all_years.groupby('time').mean([\"lat\", \"lon\"])\n",
    "# cali_seasonal_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This produces a raster for each season over time across regions\n",
    "west_coast_mean_temp_raster = west_temp_masked.resample(\n",
    "    time='QS-DEC', keep_attrs=True).mean()\n",
    "west_coast_mean_temp_raster.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "caption": "Facet plot showing seasonal mean temperature values for the several areas of interest (aoi's) including California, Washington, Oregon and Nevada and summarized over time."
   },
   "outputs": [],
   "source": [
    "# This produces a regional summary\n",
    "regional_summary = west_coast_mean_temp_raster.groupby('time').mean([\n",
    "    \"lat\", \"lon\"])\n",
    "regional_summary.plot(col=\"region\",\n",
    "                      marker=\"o\",\n",
    "                      color=\"grey\",\n",
    "                      markerfacecolor=\"purple\",\n",
    "                      markeredgecolor=\"purple\",\n",
    "                      col_wrap=2,\n",
    "                      figsize=(12, 8))\n",
    "\n",
    "plt.suptitle(\"Seasonal Temperature by Region\", y=1.03)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The data can then be easily converted to a dataframe\n",
    "regional_summary.to_dataframe()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "285.8px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}