SeaDataNet_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example Notebook for SeaDataNet Climatology data\n",
    "Author(s): [Bjorn Backeberg](mailto:backeb@gmail.com) (backeb) <br> \n",
    "Creation date: 01-Aug-2019 <br>\n",
    "Last updated:  06-Aug-2019 <br>\n",
    "\n",
    "---\n",
    "\n",
    "## Purpose\n",
    "1. Load SeaDataNet Climatology computed from the SeaDataNet V1.1 aggregated regional datasets. Data can be downloaded [here](https://www.seadatanet.org/Products#/search?from=1&to=20)\n",
    "2. Plot on map using cartopy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-06T11:39:56.745964Z",
     "start_time": "2019-08-06T11:39:56.167514Z"
    }
   },
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-06T11:39:57.170880Z",
     "start_time": "2019-08-06T11:39:57.129512Z"
    }
   },
   "outputs": [],
   "source": [
    "# open dataset for nearest user defined location\n",
    "ds1 = xr.open_dataset('data/SDN_Clim_Arctic_Temperature.nc')\n",
    "ds2 = xr.open_dataset('data/SDN_Clim_BalticSea_Temperature.nc')\n",
    "\n",
    "# display some of the metadata for ds1\n",
    "print(ds1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-06T11:39:59.122966Z",
     "start_time": "2019-08-06T11:39:58.835431Z"
    }
   },
   "outputs": [],
   "source": [
    "# load the data into variables for plotting\n",
    "lon1 = ds1.lon.values\n",
    "lat1 = ds1.lat.values\n",
    "depth1 = ds1.depth.values\n",
    "time1 = ds1.time.values\n",
    "temperature1 = ds1.Temperature.values\n",
    "\n",
    "lon2 = ds2.lon.values\n",
    "lat2 = ds2.lat.values\n",
    "depth2 = ds2.depth.values\n",
    "time2 = ds2.time.values\n",
    "temperature2 = ds2.Temperature.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot on map using cartopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-06T11:40:09.424213Z",
     "start_time": "2019-08-06T11:40:00.598870Z"
    }
   },
   "outputs": [],
   "source": [
    "tm = 6 # set time to plot\n",
    "dp = -1 # set depth to plot\n",
    "\n",
    "# instantiate the figure\n",
    "fig = plt.figure(figsize = (13, 10), dpi = 80) \n",
    "\n",
    "# set cartopy projection\n",
    "central_longitude = np.median(np.concatenate((lon1, lon2), axis = None)).round()\n",
    "central_latitude = np.median(np.concatenate((lat1, lat2), axis = None)).round()\n",
    "ax = plt.axes(projection=ccrs.NearsidePerspective(central_longitude = central_longitude,\n",
    "                                                  central_latitude = central_latitude))\n",
    "# define temperature colour axis bounds\n",
    "bounds = [np.nanmin(np.concatenate((temperature1, temperature2), axis = None)), \n",
    "          np.nanmax(np.concatenate((temperature1, temperature2), axis = None))]\n",
    "\n",
    "# plot\n",
    "cm = ax.pcolormesh(lon1, lat1, temperature1[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r, \n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "cm = ax.pcolormesh(lon2, lat2, temperature2[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r,\n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "\n",
    "plt.colorbar(cm, orientation = 'vertical').set_label('[$^\\circ$C]')\n",
    "\n",
    "# add a coastline\n",
    "coastline = cfeature.GSHHSFeature(scale = 'intermediate', edgecolor = 'none', facecolor = 'grey')\n",
    "ax.add_feature(coastline)\n",
    "\n",
    "# add a title\n",
    "plt.title(\"SeaDataNet Temperature for \"\n",
    "          +str((pd.to_datetime(time1[tm])).strftime('%B'))\n",
    "          +\", depth = \"+str(depth1[dp])+\"m\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1\n",
    "The above demo gives you an example of how to plot data from two regions on a map. We have plotted the following two data:\n",
    "* `SDN_Clim_Arctic_Temperature.nc`\n",
    "* `SDN_Clim_BalticSea_Temperature.nc`\n",
    "\n",
    "For this exercise you will need to add the `SDN_Clim_BlackSea_Temperature.nc` to the map. Your final map should look something like this: <br>\n",
    "![alt text](figs/image.png \"Logo Title Text 1\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2\n",
    "Plot the salinity from all three regions on a map. The data files you will need are:\n",
    "* `SDN_Clim_Arctic_Salinity.nc`\n",
    "* `SDN_Clim_BalticSea_Salinity.nc`\n",
    "* `SDN_Clim_BlackSea_Salinity.nc`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Solution for Exercise 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# open dataset for nearest user defined location\n",
    "ds1 = xr.open_dataset('data/SDN_Clim_Arctic_Temperature.nc')\n",
    "ds2 = xr.open_dataset('data/SDN_Clim_BalticSea_Temperature.nc')\n",
    "ds3 = xr.open_dataset('data/SDN_Clim_BlackSea_Temperature.nc')\n",
    "\n",
    "# load the data into variables for plotting\n",
    "lon1 = ds1.lon.values\n",
    "lat1 = ds1.lat.values\n",
    "depth1 = ds1.depth.values\n",
    "time1 = ds1.time.values\n",
    "temperature1 = ds1.Temperature.values\n",
    "\n",
    "lon2 = ds2.lon.values\n",
    "lat2 = ds2.lat.values\n",
    "depth2 = ds2.depth.values\n",
    "time2 = ds2.time.values\n",
    "temperature2 = ds2.Temperature.values\n",
    "\n",
    "lon3 = ds3.lon.values\n",
    "lat3 = ds3.lat.values\n",
    "depth3 = ds3.depth.values\n",
    "time3 = ds3.time.values\n",
    "temperature3 = ds3.Temperature.values\n",
    "\n",
    "tm = 6 # set time to plot\n",
    "dp = -1 # set depth to plot\n",
    "\n",
    "# instantiate the figure\n",
    "fig = plt.figure(figsize = (13, 10), dpi = 80) \n",
    "\n",
    "# set cartopy projection\n",
    "central_longitude = np.median(np.concatenate((lon1,lon2,lon3), axis = None)).round()\n",
    "central_latitude = np.median(np.concatenate((lat1,lat2,lat3), axis = None)).round()\n",
    "ax = plt.axes(projection=ccrs.NearsidePerspective(central_longitude = central_longitude,\n",
    "                                                  central_latitude = central_latitude))\n",
    "# define temperature colour axis bounds\n",
    "bounds = [np.nanmin(np.concatenate((temperature1, temperature2, temperature3), axis = None)), \n",
    "          np.nanmax(np.concatenate((temperature1, temperature2, temperature3), axis = None))]\n",
    "\n",
    "# plot\n",
    "cm = ax.pcolormesh(lon1, lat1, temperature1[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r, \n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "cm = ax.pcolormesh(lon2, lat2, temperature2[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r,\n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "cm = ax.pcolormesh(lon3, lat3, temperature3[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r,\n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "plt.colorbar(cm, orientation = 'vertical').set_label('[$^\\circ$C]')\n",
    "\n",
    "# add a coastline\n",
    "coastline = cfeature.GSHHSFeature(scale = 'intermediate', edgecolor = 'none', facecolor = 'grey')\n",
    "ax.add_feature(coastline)\n",
    "\n",
    "# add a title\n",
    "plt.title(\"SeaDataNet Temperature for \"\n",
    "          +str((pd.to_datetime(time1[tm])).strftime('%B'))\n",
    "          +\", depth = \"+str(depth1[dp])+\"m\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Solution for Exercise 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-06T12:11:35.687318Z",
     "start_time": "2019-08-06T12:11:33.935414Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# open dataset for nearest user defined location\n",
    "ds1 = xr.open_dataset('data/SDN_Clim_Arctic_Salinity.nc')\n",
    "ds2 = xr.open_dataset('data/SDN_Clim_BalticSea_Salinity.nc')\n",
    "ds3 = xr.open_dataset('data/SDN_Clim_BlackSea_Salinity.nc')\n",
    "\n",
    "# load the data into variables for plotting\n",
    "lon1 = ds1.lon.values\n",
    "lat1 = ds1.lat.values\n",
    "depth1 = ds1.depth.values\n",
    "time1 = ds1.time.values\n",
    "salinity1 = ds1.Salinity.values\n",
    "\n",
    "lon2 = ds2.lon.values\n",
    "lat2 = ds2.lat.values\n",
    "depth2 = ds2.depth.values\n",
    "time2 = ds2.time.values\n",
    "salinity2 = ds2.Salinity.values\n",
    "\n",
    "lon3 = ds3.lon.values\n",
    "lat3 = ds3.lat.values\n",
    "depth3 = ds3.depth.values\n",
    "time3 = ds3.time.values\n",
    "salinity3 = ds3.Salinity.values\n",
    "\n",
    "tm = 6 # set time to plot\n",
    "dp = -1 # set depth to plot\n",
    "\n",
    "# instantiate the figure\n",
    "fig = plt.figure(figsize = (13, 10), dpi = 80) \n",
    "\n",
    "# set cartopy projection\n",
    "central_longitude = np.median(np.concatenate((lon1,lon2,lon3), axis = None)).round()\n",
    "central_latitude = np.median(np.concatenate((lat1,lat2,lat3), axis = None)).round()\n",
    "ax = plt.axes(projection=ccrs.NearsidePerspective(central_longitude = central_longitude,\n",
    "                                                  central_latitude = central_latitude))\n",
    "# define salinity colour axis bounds\n",
    "bounds = [np.nanmin(np.concatenate((salinity1, salinity2, salinity3), axis = None)), \n",
    "          np.nanmax(np.concatenate((salinity1, salinity2, salinity3), axis = None))]\n",
    "\n",
    "# plot\n",
    "cm = ax.pcolormesh(lon1, lat1, salinity1[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r, \n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "cm = ax.pcolormesh(lon2, lat2, salinity2[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r,\n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "cm = ax.pcolormesh(lon3, lat3, salinity3[tm,dp,:,:], \n",
    "              shading = 'gourand', \n",
    "              cmap = plt.cm.Spectral_r,\n",
    "              vmin = bounds[0], vmax = bounds[1],\n",
    "              transform = ccrs.PlateCarree())\n",
    "\n",
    "plt.colorbar(cm, orientation = 'vertical').set_label('psu')\n",
    "\n",
    "# add a coastline\n",
    "coastline = cfeature.GSHHSFeature(scale = 'intermediate', edgecolor = 'none', facecolor = 'grey')\n",
    "ax.add_feature(coastline)\n",
    "\n",
    "# add a title\n",
    "plt.title(\"SeaDataNet Salinity for \"\n",
    "          +str((pd.to_datetime(time1[tm])).strftime('%B'))\n",
    "          +\", depth = \"+str(depth1[dp])+\"m\")\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.9 64-bit",
   "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.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}