{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_overlay_11a.py\nThis script illustrates the following concepts:\n    - Overlaying vectors and filled contours on a map\n    - Masking out particular areas in a map\n    - Subsetting a color map\n\nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/overlay_11.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/overlay_11_lg.png\n\nThis script shows how to overlay contours and vectors on a map,\nbut with the contours limited to specific areas, and the vectors\nnot limited.\n\nThe point of this script is to show how to mask contours against\na geographical boundary, but in a way that allows them to be drawn\nup to the boundary location. This is unlike the shapefile masking\nexamples, where grid points are set to missing if they fall\noutside a boundary, and hence you can get blocky features close\nto the boundary.\n\nWith Python's matplotlib, there are 2 general approaches to\naccomplishing this:\n\na. You can \"cover\" (i.e., \"over lay\") geographical features on\n   top of other plots using the ``zorder`` parameter to most\n   rendered objects.\n\nb. You can \"clip\" a plot object with a geographical boundary.\n\nThis example demonstrates approach (a).\n\nSee NCL_overlay_11b.py for demonstration of approach (b).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Import packages:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import xarray as xr\nimport numpy as np\n\nfrom matplotlib import pyplot as plt\n\nfrom cartopy.feature import ShapelyFeature, OCEAN, LAKES, LAND\nfrom cartopy.crs import PlateCarree\nfrom cartopy.io.shapereader import Reader as ShapeReader, natural_earth\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import cmaps as gvcmaps\nfrom geocat.viz import util as gvutil"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read in data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Open a netCDF data file using xarray default engine and load the data into xarrays, as well as extract slices for\n# ``time=0`` and the ``lev=500`` hPa level\nds = xr.open_dataset(gdf.get(\"netcdf_files/uvt.nc\")).sel(time=0, lev=500)\n\n# For convenience only, extract the U,V,T and lat and lon variables\nU = ds[\"U\"]\nV = ds[\"V\"]\nT = ds[\"T\"]\n\nlat = ds[\"lat\"]\nlon = ds[\"lon\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Construct shape boundaries:\n\nUsing Cartopy's interface to the Natural Earth Collection of shapefiles\nand geographical shape data, we construct the geographical boundaries\nthat we are interested in displaying, namely the country borders of China\nand Taiwan, the borders of Chinese provinces, and all land borders *without*\nChina or Taiwan.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Download the Natural Earth shapefile for country boundaries at 10m resolution\nshapefile = natural_earth(category='cultural',\n                          resolution='10m',\n                          name='admin_0_countries')\n\n# Sort the geometries in the shapefile into Chinese/Taiwanese or other\ncountry_geos = []\nother_land_geos = []\nfor record in ShapeReader(shapefile).records():\n    if record.attributes['ADMIN'] in ['China', 'Taiwan']:\n        country_geos.append(record.geometry)\n    else:\n        other_land_geos.append(record.geometry)\n\n# Define map projection to allow Cartopy to transform ``lat`` and ``lon`` values accurately into points on the\n# matplotlib plot canvas.\nprojection = PlateCarree()\n\n# Define a Cartopy Feature for the country borders and the land mask (i.e.,\n# all other land) from the shapefile geometries, so they can be easily plotted\ncountries = ShapelyFeature(country_geos,\n                           crs=projection,\n                           facecolor='none',\n                           edgecolor='black',\n                           lw=1.5)\nland_mask = ShapelyFeature(other_land_geos,\n                           crs=projection,\n                           facecolor='white',\n                           edgecolor='none')\n\n# Download the Natural Earth shapefile for the states/provinces at 10m resolution\nshapefile = natural_earth(category='cultural',\n                          resolution='10m',\n                          name='admin_1_states_provinces')\n\n# Extract the Chinese province borders\nprovince_geos = [\n    record.geometry\n    for record in ShapeReader(shapefile).records()\n    if record.attributes['admin'] == 'China'\n]\n\n# Define a Cartopy Feature for the province borders, so they can be easily plotted\nprovinces = ShapelyFeature(province_geos,\n                           crs=projection,\n                           facecolor='none',\n                           edgecolor='black',\n                           lw=0.25)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Generate figure (set its size (width, height) in inches) and axes using Cartopy\nfig = plt.figure(figsize=(10, 10))\nax = plt.axes(projection=projection)\n\nax.set_extent([100, 145, 15, 55], crs=projection)\n\n# Define the contour levels\nclevs = np.arange(228, 273, 4, dtype=float)\n\n# Import an NCL colormap, truncating it by using geocat.viz.util convenience function\nnewcmp = gvutil.truncate_colormap(gvcmaps.BkBlAqGrYeOrReViWh200,\n                                  minval=0.1,\n                                  maxval=0.6,\n                                  n=len(clevs))\n\n# Draw the temperature contour plot with the subselected colormap\n# (Place the zorder of the contour plot at the lowest level)\ncf = ax.contourf(lon, lat, T, levels=clevs, cmap=newcmp, zorder=1)\n\n# Draw horizontal color bar\ncax = plt.axes((0.14, 0.08, 0.74, 0.02))\ncbar = plt.colorbar(cf,\n                    ax=ax,\n                    cax=cax,\n                    ticks=clevs[1:-1],\n                    drawedges=True,\n                    orientation='horizontal')\ncbar.ax.tick_params(labelsize=12)\n\n# Add the land mask feature on top of the contour plot (higher zorder)\nax.add_feature(land_mask, zorder=2)\n\n# Add the OCEAN and LAKES features on top of the contour plot\nax.add_feature(OCEAN.with_scale('50m'), edgecolor='black', lw=1, zorder=2)\nax.add_feature(LAKES.with_scale('50m'), edgecolor='black', lw=1, zorder=2)\n\n# Add the country and province features (which are transparent) on top\nax.add_feature(countries, zorder=3)\nax.add_feature(provinces, zorder=3)\n\n# Draw the wind quiver plot on top of everything else\nQ = ax.quiver(lon,\n              lat,\n              U,\n              V,\n              color='black',\n              width=.003,\n              scale=600.,\n              headwidth=3.75,\n              zorder=4)\n\n# Draw the key for the quiver plot\nrect = plt.Rectangle((142, 52),\n                     3,\n                     3,\n                     facecolor='mediumorchid',\n                     edgecolor=None,\n                     zorder=4)\nax.add_patch(rect)\nax.quiverkey(Q,\n             0.9675,\n             0.95,\n             30,\n             '30',\n             labelpos='N',\n             color='black',\n             coordinates='axes',\n             fontproperties={'size': 14},\n             labelsep=0.1)\n\n# Add a text box to indicate the pressure level\nprops = dict(facecolor='white', edgecolor='none', alpha=0.8)\nax.text(105,\n        52.7,\n        '500hPa',\n        transform=projection,\n        fontsize=18,\n        ha='center',\n        va='center',\n        color='mediumorchid',\n        bbox=props)\n\n# Use geocat.viz.util convenience function to set axes tick values\ngvutil.set_axes_limits_and_ticks(ax,\n                                 xticks=[100, 120, 140],\n                                 yticks=[20, 30, 40, 50])\n\n# Use geocat.viz.util convenience function to make plots look like NCL plots by using latitude, longitude tick labels\ngvutil.add_lat_lon_ticklabels(ax)\n\n# Use geocat.viz.util convenience function to add minor and major tick lines\ngvutil.add_major_minor_ticks(ax,\n                             x_minor_per_major=4,\n                             y_minor_per_major=5,\n                             labelsize=18)\n\n# Use geocat.viz.util convenience function to add main title as well as titles to left and right of the plot axes.\ngvutil.set_titles_and_labels(ax,\n                             lefttitle=\"Temp\",\n                             lefttitlefontsize=20,\n                             righttitle=\"Wind\",\n                             righttitlefontsize=20)\n\n# Show the plot\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "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": 0
}