netCDF, xarray, and cartopy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "z3hprDvHs9jN"
   },
   "source": [
    "# Geospatial data analysis and visualization using xarray, matplotlib, cartopy, and MetPy\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VJZ638JUs9jO"
   },
   "source": [
    "So far, we have used `matplotlib` to visualize data.  In this exercise, we will put it all together: \n",
    "- read files\n",
    "- subset files using fancy indexing\n",
    "- plot the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EaQR3JP3s9jP"
   },
   "source": [
    "We're going to use a new python package called `xarray`.  It is a high level interface for netCDF file reading and writing, that will make your life easier.  Let's go!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "VyrPkvAls9jQ"
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hIS-tugxs9jW"
   },
   "source": [
    "Ahhhhhhhh....all better.  Let's get to work!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iie9lod7s9jW"
   },
   "source": [
    "## Using xarray to read a netCDF file\n",
    "\n",
    "netCDF4 is a common dataset for storing gridded binary data in atmospheric sciences.  It is a self-describing data format, meaning that it contains data and all of the coordinates necessary to use the data.  We will start with a simple example - maps of surface temperature anomalies from NASA GISS (from a server at NOAA ESRL).  This particular dataset is online and available for streaming through a service called OPENDaP - which means we don't even have to download the data.\n",
    "\n",
    "It can't be much easier to read data than this...xarray handles a lot of the dirty work for you.  We can load both local files, as well as files on the internet like this OPENDaP file.  Xarray allows you to either give the local file path, or the web site!\n",
    "\n",
    "Why xarray? It allows interfaces to self-describing data (such as netCDF) that is structured so you know exactly what you are getting, provides tools for easy analysis of the data.\n",
    "\n",
    "<img src=\"https://github.com/pydata/xarray/raw/master/doc/_static/dataset-diagram.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "owEkZSYas9jX"
   },
   "outputs": [],
   "source": [
    "nc=xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/gistemp/combined/250km/air.2x2.250.mon.anom.comb.nc')\n",
    "nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CPPDhg8Ws9jZ"
   },
   "outputs": [],
   "source": [
    "ncvar = nc['air']\n",
    "ncvar\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "paD9gBtnvDg2"
   },
   "outputs": [],
   "source": [
    "air_tavg = ncvar.mean(dim='time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "O820Mb-IwMng"
   },
   "outputs": [],
   "source": [
    "air_tavg.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MGsZXwEXs9ja"
   },
   "source": [
    "This data is a gridded time series at 2 degree resolution of monthly surface temperature anomalies starting in 1880.  Let's average over all space dimensions (lat - axis 1, and lon - axis 2).  We can use `np.mean` and its `axis` keyword (very handy) for this purpose.\n",
    "\n",
    "This will yield a time series of globally averaged temperature!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1Apy21Gxs9jb"
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(ncvar.time,np.mean(ncvar,axis=(1,2)))\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Temperature anomaly deg C')\n",
    "plt.title(\"It's getting hot up in here!\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fDS1bT7xxQ64"
   },
   "outputs": [],
   "source": [
    "ncvar.mean(dim=['lat','lon']).plot()\n",
    "plt.xlabel('Years since some time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "z1Ky1efjs9jd"
   },
   "source": [
    "We can use the select tool to get a subset in a box (find closest index values of lon and lat) so that we can subset the data and grab the closest point to Champaign-Urbana.  We can give it a list of points.  Here we will give it one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "eS0GQCB7s9je"
   },
   "outputs": [],
   "source": [
    "nc_cmi=nc.sel(lon=-88.9+360.,lat=40., method='nearest')\n",
    "nc_cmi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hDV2OBK8s9jh"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "var='air'\n",
    "\n",
    "plt.plot(nc_cmi.time,nc_cmi[var])\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Temperature anomaly deg C')\n",
    "plt.title(\"It's not getting as hot in Champaign\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QJlM4f0Gyx8D"
   },
   "outputs": [],
   "source": [
    "nc['air'].sel(lon=-88.9+360.,lat=40., method='nearest').plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "g6rD4ldPs9jj"
   },
   "source": [
    "## Saving to file - easy as np.pi()\n",
    "Want to save the file as a netCDF file?  No problem!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "UhnTLCuYs9jk"
   },
   "outputs": [],
   "source": [
    "nc_cmi.to_netcdf('nc_cmi.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZmN6TUr4s9jl"
   },
   "source": [
    "## Calculating time averages\n",
    "\n",
    "Let's say we want to average the monthly time series data into annually averaged data.\n",
    "\n",
    "There are a number of ways to do this.  `xarray` offers time sampling capabilities, similar to `pandas`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "blV1n2-ws9jm"
   },
   "outputs": [],
   "source": [
    "nc_cmi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NV2EIo3Vs9jx"
   },
   "source": [
    "`xarray` and `pandas` share the same interface to resample and group time series conveniently.  The documentation is available at: http://xarray.pydata.org/en/stable/time-series.html#resampling-and-grouped-operations.  The codes for resampling are the same as `pandas`.  See http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cKo5qHn6s9jy"
   },
   "outputs": [],
   "source": [
    "nc_cmi3=nc_cmi.resample(time='AS').mean()\n",
    "nc_cmi3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sMV2Wxh_s9j1"
   },
   "source": [
    "How easy is that?  np.pi()?  You can also resample with other time frequencies, or in space, or change how you do the calculation (i.e., calculate the median instead of the mean).\n",
    "\n",
    "Now save to a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "99y2t9Ous9j2"
   },
   "outputs": [],
   "source": [
    "site='Champaign_annualavgs'\n",
    "\n",
    "nc_cmi3.to_netcdf(site+'_data2.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "J5XtHQkls9j4"
   },
   "source": [
    "Make a (nice) plot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J_aNd2NVs9j4"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(11,8.5)) #create a new figure\n",
    "\n",
    "plt.plot(nc_cmi['time'],np.squeeze(nc_cmi['air']),'b',alpha=0.5)\n",
    "plt.plot(nc_cmi3['time'],np.squeeze(nc_cmi3['air']),'r',linewidth=2.0)\n",
    "plt.legend(['Monthly averages','Annual averages'])\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Temperature Anomaly (degrees C)')\n",
    "plt.title('GISTEMP Temperature Anomalies near Champaign, IL')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rFPpGRa1wRwh"
   },
   "source": [
    "## Collapsing dimensions to make a Hovmoller Diagram\n",
    "\n",
    "It's also easy to slice and dice the data to make a Hovmoller Diagram (time,longitude plot averaging over latitude)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rFPpGRa1wRwh"
   },
   "outputs": [],
   "source": [
    "temp1= nc['air'].resample(time='AS').mean()\n",
    "temp1.mean(dim='lat').plot(vmin=-4,vmax=4,cmap='seismic')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GDC7MZJUziWx"
   },
   "source": [
    "Mapping with `cartopy`\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "eEbFNbIToDli"
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "\n",
    "import xarray as xr\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tw1mJr5bziW2"
   },
   "source": [
    "## Introduction to mapping with `cartopy`.\n",
    "\n",
    "Python's `cartopy` is a package that is designed to enable `matplotlib` handle geographic calculations and visualize maps with geographic boundaries.  We can then use `matplotlib` to make maps, and off we go!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IlAFQTSYziW6"
   },
   "source": [
    "Cartopy has exposed an interface to enable easy map creation using matplotlib. Creating a basic map is as simple as telling matplotlib to use a specific map projection, and then adding some coastlines to the axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "XzYxubE7ziW7"
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.add_feature(cfeature.COASTLINE)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Ob9olZhoziW9"
   },
   "source": [
    "A list of the available projections to be used with `matplotlib` can be found on the Cartopy projection list page.\n",
    "\n",
    "The line `plt.axes(projection=ccrs.PlateCarree())` sets up a `GeoAxes` instance which exposes a variety of other map related methods, in the case of the previous example, we used the `coastlines()` method to add coastlines to the map.\n",
    "\n",
    "Lets create another map in a different projection, and make use of the `stock_img()` method to add an underlay image to the map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KtCWKwpTziW-"
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(11,11))\n",
    "ax = plt.axes(projection=ccrs.Mollweide())\n",
    "ax.coastlines()\n",
    "\n",
    "ax.stock_img()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mu3PQwQrziXA"
   },
   "source": [
    "## Adding data to the map\n",
    "Once you have the map just the way you want it, data can be added to it in exactly the same way as with normal matplotlib axes. By default, the coordinate system of any data added to a GeoAxes is the same as the coordinate system of the GeoAxes itself, to control which coordinate system that the given data is in, you can add the `transform` keyword with an appropriate `cartopy.crs.CRS` instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tuMcDjBCziXB"
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "ax = plt.axes(projection=ccrs.Orthographic(central_longitude=20,central_latitude=40))\n",
    "ax.coastlines()\n",
    "\n",
    "ny_lon, ny_lat = -75, 43\n",
    "delhi_lon, delhi_lat = 77.23, 28.61\n",
    "\n",
    "plt.plot([ny_lon, delhi_lon], [ny_lat, delhi_lat],\n",
    "         color='blue', linewidth=2, marker='o',\n",
    "         transform=ccrs.Geodetic(),\n",
    "         )\n",
    "\n",
    "plt.plot([ny_lon, delhi_lon], [ny_lat, delhi_lat],\n",
    "         color='gray', linestyle='--',\n",
    "         transform=ccrs.PlateCarree(),\n",
    "         )\n",
    "\n",
    "plt.text(ny_lon - 3, ny_lat - 12, 'New York',\n",
    "         horizontalalignment='right',\n",
    "         transform=ccrs.Geodetic())\n",
    "\n",
    "plt.text(delhi_lon + 3, delhi_lat - 12, 'Delhi',\n",
    "         horizontalalignment='left',\n",
    "         transform=ccrs.Geodetic())\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "J0nrzJ6CziXD"
   },
   "source": [
    "# More advanced mapping with `cartopy` and `matplotlib`\n",
    "\n",
    "From the outset, `cartopy`’s purpose has been to simplify and improve the quality of mapping visualisations available for scientific data. Thanks to the simplicity of the cartopy interface, in many cases the hardest part of producing such visualisations is getting hold of the data in the first place. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aUy1tC_JziXE"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from cartopy import config\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "\n",
    "# get the path of the file. \n",
    "sst_data = xr.open_dataset('https://www.ncei.noaa.gov/thredds/dodsC/ncFC/fc-oisst-daily-avhrr-only-dly/OISST_Daily_AVHRR-only_Feature_Collection_best.ncd').sel(time='2020-08-01').isel(zlev=0)\n",
    "\n",
    "sst_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6dCFPeCeziXH"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "\n",
    "plt.pcolormesh(sst_data['lon'], sst_data['lat'], sst_data['sst'].isel(time=0),\n",
    "             transform=ccrs.PlateCarree(), cmap='viridis')\n",
    "\n",
    "plt.title(sst_data.title+'\\n'+str(sst_data['time'].values))\n",
    "ax.coastlines()\n",
    "\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gZBQHtEBziXK"
   },
   "source": [
    "## Vector plotting\n",
    "\n",
    "Cartopy comes with powerful vector field plotting functionality. There are 3 distinct options for visualising vector fields: quivers, barbs, and streamplots, each with their own benefits for displaying certain vector field forms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-K_VysBZziXL"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import cartopy\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "\n",
    "def sample_data(shape=(20, 30)):\n",
    "    \"\"\"\n",
    "    Returns ``(x, y, u, v, crs)`` of some vector data\n",
    "    computed mathematically. The returned crs will be a rotated\n",
    "    pole CRS, meaning that the vectors will be unevenly spaced in\n",
    "    regular PlateCarree space.\n",
    "\n",
    "    \"\"\"\n",
    "    crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)\n",
    "\n",
    "    x = np.linspace(311.9, 391.1, shape[1])\n",
    "    y = np.linspace(-23.6, 24.8, shape[0])\n",
    "\n",
    "    x2d, y2d = np.meshgrid(x, y)\n",
    "    u = 10 * (2 * np.cos(2 * np.deg2rad(x2d) + 3 * np.deg2rad(y2d + 30)) ** 2)\n",
    "    v = 20 * np.cos(6 * np.deg2rad(x2d))\n",
    "\n",
    "    return x, y, u, v, crs\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "\n",
    "ax = plt.axes(projection=ccrs.Orthographic(-10, 45))\n",
    "\n",
    "ax.add_feature(cartopy.feature.OCEAN, zorder=0)\n",
    "ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')\n",
    "\n",
    "ax.set_global()\n",
    "ax.gridlines()\n",
    "\n",
    "x, y, u, v, vector_crs = sample_data()\n",
    "ax.streamplot(x, y, u, v, transform=vector_crs)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XSOgPaLrziXO"
   },
   "source": [
    "## Putting it back together\n",
    "\n",
    "Now, let's work with our global temperature anomaly dataset.  We will use xarray to select some times, and visualize the anomalies on a two panel plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6RYXdiTiziXP"
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import matplotlib.ticker as mticker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Iv5Y8wADziXS"
   },
   "outputs": [],
   "source": [
    "nc=xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/gistemp/combined/250km/air.2x2.250.mon.anom.comb.nc')\n",
    "nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QwoUT435ziXV"
   },
   "outputs": [],
   "source": [
    "#Select times, make lats and lons 2D so basemap knows where to plot each point\n",
    "\n",
    "jan1976 = nc.sel(time=pd.datetime(1976,1,1))\n",
    "jan2018 = nc.sel(time=pd.datetime(2018,1,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VZmENg0qziXX"
   },
   "source": [
    "Here, let's plot the anomalies from these two time periods on a rectangular map projection called \"PlateCarree\" in a two panel plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6MlLoi4PCLHa"
   },
   "outputs": [],
   "source": [
    "jan2018['air'].plot.pcolormesh()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QjOVNoX4ziXY"
   },
   "outputs": [],
   "source": [
    "# select data we want\n",
    "data = jan1976['air']\n",
    "data2 = jan2018['air']\n",
    "\n",
    "# create figure, axes instances.\n",
    "fig = plt.figure(figsize=(12, 3))\n",
    "\n",
    "ax = plt.subplot(1, 2, 1, projection=ccrs.PlateCarree(central_longitude=180.))\n",
    "data.plot.contourf(ax=ax, transform=ccrs.PlateCarree(), levels=np.linspace(-12, 12, 25));\n",
    "ax.set_global(); \n",
    "ax.coastlines();\n",
    "#set up the gridlines\n",
    "gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,\n",
    "                  linewidth=1, color='gray', alpha=0.5, linestyle='-')\n",
    "#set where the gridlines go\n",
    "gl.xlocator = mticker.FixedLocator(np.arange(-180,181,15))\n",
    "gl.ylocator = mticker.FixedLocator(np.arange(-90,91,15))\n",
    "\n",
    "ax2 = plt.subplot(1, 2, 2, projection=ccrs.PlateCarree())\n",
    "data2.plot.contourf(ax=ax2, transform=ccrs.PlateCarree(), levels=np.linspace(-12, 12, 25));\n",
    "ax2.set_global(); \n",
    "ax2.coastlines();\n",
    "\n",
    "gl2 = ax2.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,\n",
    "                  linewidth=1, color='gray', alpha=0.5, linestyle='-')\n",
    "\n",
    "gl2.ylocator = mticker.FixedLocator(np.arange(-90,91,15))\n",
    "gl2.xlocator = mticker.FixedLocator(np.arange(-180,181,15))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VQK2T9J_ziXa"
   },
   "source": [
    "Let's zoom in to the Eastern US - and add some fancy stuff like country and state borders, and lakes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JRubn-0VziXb"
   },
   "outputs": [],
   "source": [
    "import cartopy.feature as cfeature\n",
    "\n",
    "# select data we want\n",
    "data = jan1976['air']\n",
    "data2 = jan2018['air']\n",
    "\n",
    "# create figure, axes instances.\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "\n",
    "cartopy_proj = projection=ccrs.Mercator()\n",
    "\n",
    "ax = plt.subplot(1, 2, 1, projection=cartopy_proj)\n",
    "data.plot.contourf(ax=ax, transform=ccrs.PlateCarree(), levels=np.linspace(-12, 12, 25));\n",
    "ax.set_extent([-100,-80,30.,50.], crs=ccrs.PlateCarree()) \n",
    "ax.coastlines()\n",
    "# Create a feature for States/Admin 1 regions at 1:50m from Natural Earth\n",
    "states_provinces = cfeature.NaturalEarthFeature(\n",
    "    category='cultural',\n",
    "    name='admin_1_states_provinces_lines',\n",
    "    scale='50m',\n",
    "    facecolor='none')\n",
    "ax.add_feature(states_provinces, edgecolor='gray')\n",
    "ax.add_feature(cfeature.LAKES, edgecolor='gray')\n",
    "ax.add_feature(cfeature.BORDERS, edgecolor='gray')\n",
    "#set up the gridlines\n",
    "gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,\n",
    "                  linewidth=2, color='gray', alpha=0.5, linestyle='-')\n",
    "#set where the gridlines go\n",
    "gl.xlocator = mticker.FixedLocator(np.arange(-180,181,5))\n",
    "gl.ylocator = mticker.FixedLocator(np.arange(-90,90,5))\n",
    "\n",
    "ax2 = plt.subplot(1, 2, 2, projection=cartopy_proj)\n",
    "data2.plot.contourf(ax=ax2, transform=ccrs.PlateCarree(), levels=np.linspace(-12, 12, 25));\n",
    "ax2.set_extent([-100,-80,30.,50.], crs=ccrs.PlateCarree()) \n",
    "ax2.coastlines()\n",
    "# Create a feature for States/Admin 1 regions at 1:50m from Natural Earth\n",
    "states_provinces = cfeature.NaturalEarthFeature(\n",
    "    category='cultural',\n",
    "    name='admin_1_states_provinces_lines',\n",
    "    scale='50m',\n",
    "    facecolor='none')\n",
    "ax2.add_feature(states_provinces, edgecolor='gray')\n",
    "ax2.add_feature(cfeature.LAKES, edgecolor='gray')\n",
    "ax2.add_feature(cfeature.BORDERS, edgecolor='gray')\n",
    "gl2 = ax2.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,\n",
    "                  linewidth=2, color='gray', alpha=0.5, linestyle='-')\n",
    "\n",
    "gl2.ylocator = mticker.FixedLocator(np.arange(-90,90,5))\n",
    "gl2.xlocator = mticker.FixedLocator(np.arange(-180,181,5))\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "g6rD4ldPs9jj",
    "ZmN6TUr4s9jl",
    "GDC7MZJUziWx",
    "tw1mJr5bziW2",
    "mu3PQwQrziXA",
    "gZBQHtEBziXK"
   ],
   "name": "netCDF, xarray, and cartopy.ipynb",
   "provenance": [
    {
     "file_id": "12MtgIQgUDHT3uYxyBb6O_y6HuWEsyalB",
     "timestamp": 1601071719817
    },
    {
     "file_id": "1MQ7m8k_MR1pu5Gx1Qn6yLsTPd01ccu11",
     "timestamp": 1552409164918
    }
   ]
  },
  "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}