{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_sat_1.py\nThis script illustrates the following concepts:\n   - Creating an orthographic projection\n   - Drawing line contours over a satellite map\n   - Manually labeling contours\n   - Transforming coordinates\n\nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/sat_1.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/sat_1_lg.png\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 matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\nimport cartopy.feature as cfeature\nimport numpy as np\nfrom sklearn.cluster import DBSCAN\nimport warnings\n\nimport geocat.datafiles as gdf\nimport geocat.viz.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\n# load the data into xarrays\nds = xr.open_dataset(gdf.get(\"netcdf_files/slp.1963.nc\"), decode_times=False)\n\n# Get data from the 24th timestep\npressure = ds.slp[24, :, :]\n\n# Translate short values to float values\npressure = pressure.astype('float64')\n\n# Convert Pa to hPa data\npressure = pressure * 0.01\n\n# Fix the artifact of not-shown-data around 0 and 360-degree longitudes\nwrap_pressure = gvutil.xr_add_cyclic_longitudes(pressure, \"lon\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define a helper function to find local extrema\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def findLocalExtrema(da, highVal=0, lowVal=1000, eType='Low'):\n    \"\"\"\n    Utility function to find local low/high field variable coordinates on a contour map. To classify as a local high, the data\n    point must be greater than highVal, and to classify as a local low, the data point must be less than lowVal.\n\n    Args:\n        da: (:class:`xarray.DataArray`):\n            Xarray data array containing the lat, lon, and field variable (ex. pressure) data values\n        highVal (:class:`int`):\n            Data value that the local high must be greater than to qualify as a \"local high\" location.\n            Default highVal is 0.\n        lowVal (:class:`int`):\n            Data value that the local low must be less than to qualify as a \"local low\" location.\n            Default lowVal is 1000.\n        eType (:class:`str`):\n            'Low' or 'High'\n            Determines which extrema are being found- minimum or maximum, respectively.\n            Default eType is 'Low'.\n    Returns:\n        clusterExtremas (:class:`list`):\n            List of coordinate tuples in GPS form (lon in degrees, lat in degrees)\n            that specify local low/high locations\n    \"\"\"\n\n    # Create a 2D array of coordinates in the same shape as the field variable data\n    # so each coordinate is easily mappable to a data value\n    # ex:\n    # (1, 1), (2, 1), (3, 1)\n    # (1, 2)................\n    # (1, 3)................\n    lons, lats = np.meshgrid(np.array(da.lon), np.array(da.lat))\n    coordarr = np.dstack((lons, lats))\n\n    # Find all zeroes that also qualify as low or high values\n    extremacoords = []\n\n    if eType == 'Low':\n        coordlist = np.argwhere(da.data < lowVal)\n        extremacoords = [tuple(coordarr[x[0]][x[1]]) for x in coordlist]\n    if eType == 'High':\n        coordlist = np.argwhere(da.data < lowVal)\n        extremacoords = [tuple(coordarr[x[0]][x[1]]) for x in coordlist]\n\n    if extremacoords == []:\n        if eType == 'Low':\n            warnings.warn(\n                'No local extrema with data value less than given lowVal')\n            return []\n        if eType == 'High':\n            warnings.warn(\n                'No local extrema with data value greater than given highVal')\n            return []\n\n    # Clean up noisy data to find actual extrema\n\n    # Use Density-based spatial clustering of applications with noise\n    # to cluster and label coordinates\n    db = DBSCAN(eps=10, min_samples=1)\n    new = db.fit(extremacoords)\n    labels = new.labels_\n\n    # Create an dictionary of values with key being coordinate\n    # and value being cluster label.\n    coordsAndLabels = {label: [] for label in labels}\n    for label, coord in zip(labels, extremacoords):\n        coordsAndLabels[label].append(coord)\n\n    # Initialize array of coordinates to be returned\n    clusterExtremas = []\n\n    # Iterate through the coordinates in each cluster\n    for key in coordsAndLabels:\n\n        # Create array to hold all the field variable values for that cluster\n        datavals = []\n        for coord in coordsAndLabels[key]:\n\n            # Find pressure data at that coordinate\n            cond = np.logical_and(coordarr[:, :, 0] == coord[0],\n                                  coordarr[:, :, 1] == coord[1])\n            x, y = np.where(cond)\n            datavals.append(da.data[x[0]][y[0]])\n\n        # Find the index of the smallest/greatest field variable value of each cluster\n        if eType == 'Low':\n            index = np.argmin(np.array(datavals))\n        if eType == 'High':\n            index = np.argmax(np.array(datavals))\n\n        # Append the coordinate corresponding to that index to the array to be returned\n        clusterExtremas.append(\n            (coordsAndLabels[key][index][0], coordsAndLabels[key][index][1]))\n\n    return clusterExtremas"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define a helper function that will plot contour labels\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plotCLabels(da,\n                contours,\n                transform,\n                ax,\n                proj,\n                clabel_locations=[],\n                fontsize=12,\n                whitebbox=False,\n                horizontal=False):\n    \"\"\"\n    Utility function to plot contour labels by passing in a coordinate to the clabel function.\n    This allows the user to specify the exact locations of the labels, rather than having matplotlib\n    plot them automatically.\n\n    Args:\n        da: (:class:`xarray.DataArray`):\n            Xarray data array containing the lat, lon, and field variable data values.\n        contours (:class:`cartopy.mpl.contour.GeoContourSet`):\n            Contour set that is being labeled.\n        transform (:class:`cartopy._crs`):\n            Instance of CRS that represents the source coordinate system of coordinates.\n            (ex. ccrs.Geodetic()).\n        ax (:class:`matplotlib.pyplot.axis`):\n            Axis containing the contour set.\n        proj (:class:`cartopy.crs`):\n            Projection 'ax' is defined by.\n            This is the instaance of CRS that the coordinates will be transformed to.\n        clabel_locations (:class:`list`):\n            List of coordinate tuples in GPS form (lon in degrees, lat in degrees)\n            that specify where the contours with regular field variable values should be plotted.\n        fontsize (:class:`int`):\n            Font size of contour labels.\n        whitebbox (:class:`bool`):\n            Setting this to \"True\" will cause all labels to be plotted with white backgrounds\n        horizontal (:class:`bool`):\n            Setting this to \"True\" will cause the contour labels to be horizontal.\n    Returns:\n        cLabels (:class:`list`):\n            List of text instances of all contour labels\n    \"\"\"\n\n    # Initialize empty array that will be filled with contour label text objects and returned\n    cLabels = []\n\n    # Plot any regular contour levels\n    if clabel_locations != []:\n        clevelpoints = proj.transform_points(\n            transform, np.array([x[0] for x in clabel_locations]),\n            np.array([x[1] for x in clabel_locations]))\n        transformed_locations = [(x[0], x[1]) for x in clevelpoints]\n        ax.clabel(contours,\n                  manual=transformed_locations,\n                  inline=True,\n                  fontsize=fontsize,\n                  colors='black',\n                  fmt=\"%.0f\")\n        [cLabels.append(txt) for txt in contours.labelTexts]\n\n        if horizontal is True:\n            [txt.set_rotation('horizontal') for txt in contours.labelTexts]\n\n    if whitebbox is True:\n        [\n            txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=2))\n            for txt in cLabels\n        ]\n\n    return cLabels"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define a helper function that will plot contour labels\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plotELabels(da,\n                contours,\n                transform,\n                ax,\n                proj,\n                clabel_locations=[],\n                eType='Low',\n                fontsize=22,\n                horizontal=True,\n                whitebbox=False):\n    \"\"\"\n    Utility function to plot contour labels. High/Low contour labels will be plotted using text boxes for more accurate label values\n    and placement.\n\n    Args:\n        da: (:class:`xarray.DataArray`):\n            Xarray data array containing the lat, lon, and field variable data values.\n        contours (:class:`cartopy.mpl.contour.GeoContourSet`):\n            Contour set that is being labeled.\n        transform (:class:`cartopy._crs`):\n            Instance of CRS that represents the source coordinate system of coordinates.\n            (ex. ccrs.Geodetic()).\n        ax (:class:`matplotlib.pyplot.axis`):\n            Axis containing the contour set.\n        proj (:class:`cartopy.crs`):\n            Projection 'ax' is defined by.\n            This is the instaance of CRS that the coordinates will be transformed to.\n        clabel_locations (:class:`list`):\n            List of coordinate tuples in GPS form (lon in degrees, lat in degrees)\n            that specify where the contour labels should be plotted.\n        type (:class:`list`):\n            'high' or 'low'\n            High contour labels will be plotted with an H\n            Low contour labels will be plotted with an L\n        fontsize (:class:`int`):\n            Font size of regular contour labels.\n        horizontal (:class:`bool`):\n            Setting this to \"True\" will cause the contour labels to be horizontal.\n        whitebbox (:class:`bool`):\n            Setting this to \"True\" will cause all labels to be plotted with white backgrounds\n    Returns:\n        extremaLabels (:class:`list`):\n            List of text instances of all contour labels\n    \"\"\"\n\n    # Create array of coordinates in the same shape as field variable data\n    # so each coordinate can be easily mapped to its data value.\n    # ex:\n    # (1, 1), (2, 1), (3, 1)\n    # (1, 2)................\n    # (1, 3)................\n    lons, lats = np.meshgrid(np.array(da.lon), np.array(da.lat))\n    coordarr = np.dstack((lons, lats))\n\n    # Initialize empty array that will be filled with contour label text objects and returned\n    extremaLabels = []\n\n    # Plot any low contour levels\n    clabel_points = proj.transform_points(\n        transform, np.array([x[0] for x in clabel_locations]),\n        np.array([x[1] for x in clabel_locations]))\n    transformed_locations = [(x[0], x[1]) for x in clabel_points]\n\n    for x in range(len(transformed_locations)):\n\n        try:\n            # Find field variable data at that coordinate\n            coord = clabel_locations[x]\n            cond = np.logical_and(coordarr[:, :, 0] == coord[0],\n                                  coordarr[:, :, 1] == coord[1])\n            z, y = np.where(cond)\n            p = int(round(da.data[z[0]][y[0]]))\n\n            if eType == 'High':\n                lab = plt.text(transformed_locations[x][0],\n                               transformed_locations[x][1],\n                               \"H$_{\" + str(p) + \"}$\",\n                               fontsize=fontsize,\n                               horizontalalignment='center',\n                               verticalalignment='center')\n            elif eType == 'Low':\n                lab = plt.text(transformed_locations[x][0],\n                               transformed_locations[x][1],\n                               \"L$_{\" + str(p) + \"}$\",\n                               fontsize=fontsize,\n                               horizontalalignment='center',\n                               verticalalignment='center')\n\n            if horizontal is True:\n                lab.set_rotation('horizontal')\n\n            extremaLabels.append(lab)\n\n        except:\n            continue\n\n    if whitebbox is True:\n        [\n            txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=2))\n            for txt in extremaLabels\n        ]\n\n    return extremaLabels"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set figure size\nfig = plt.figure(figsize=(8, 8))\n\n# Set global axes with an orthographic projection\nproj = ccrs.Orthographic(central_longitude=270, central_latitude=45)\nax = plt.axes(projection=proj)\nax.set_global()\n\n# Add land, coastlines, and ocean features\nax.add_feature(cfeature.LAND, facecolor='lightgray')\nax.add_feature(cfeature.COASTLINE, linewidth=.5)\nax.add_feature(cfeature.OCEAN, facecolor='lightcyan')\nax.add_feature(cfeature.BORDERS, linewidth=.5)\nax.add_feature(cfeature.LAKES,\n               facecolor='lightcyan',\n               edgecolor='black',\n               linewidth=.5)\n\n# Make array of the contour levels that will be plotted\ncontours = np.arange(948, 1072, 4)\n\n# Plot contour data\np = wrap_pressure.plot.contour(ax=ax,\n                               transform=ccrs.PlateCarree(),\n                               linewidths=0.5,\n                               levels=contours,\n                               cmap='black',\n                               add_labels=False)\n\n# regular pressure contour levels- These values were found by setting\n# 'manual' argument in ax.clabel call to 'True' and then hovering mouse\n# over desired location of countour label to find coordinate\n# (which can be found in bottom left of figure window).\nregularCLabels = [(176.4, 34.63), (-150.46, 42.44), (-142.16, 28.5),\n                  (-134.12, 16.32), (-108.9, 17.08), (-98.17, 15.6),\n                  (-108.73, 42.19), (-111.25, 49.66), (-127.83, 41.93),\n                  (-92.49, 25.64), (-77.29, 29.08), (-77.04, 16.42),\n                  (-95.93, 57.59), (-156.05, 84.47), (-17.83, 82.52),\n                  (-76.3, 41.99), (-48.89, 41.45), (-33.43, 37.55),\n                  (-46.98, 17.17), (1.79, 63.67), (-58.78, 67.05),\n                  (-44.78, 53.68), (-69.69, 53.71), (-78.02, 52.22),\n                  (-16.91, 44.33), (-95.72, 35.17), (-102.69, 73.62)]\n\n# low pressure contour levels- these will be plotted\n# as a subscript to an 'L' symbol.\nlowCLabels = findLocalExtrema(pressure, eType='Low', highVal=1040, lowVal=975)\n\n# Plot Clabels\nplotCLabels(pressure,\n            p,\n            ccrs.Geodetic(),\n            ax,\n            proj,\n            clabel_locations=regularCLabels)\nplotELabels(pressure,\n            p,\n            ccrs.Geodetic(),\n            ax,\n            proj,\n            clabel_locations=lowCLabels,\n            eType='Low')\n\n# Use gvutil function to set title and subtitles\ngvutil.set_titles_and_labels(ax,\n                             maintitle=r\"$\\bf{SLP}$\" + \" \" + r\"$\\bf{1963,}$\" +\n                             \" \" + r\"$\\bf{January}$\" + \" \" + r\"$\\bf{24th}$\",\n                             maintitlefontsize=20,\n                             lefttitle=\"mean Daily Sea Level Pressure\",\n                             lefttitlefontsize=16,\n                             righttitle=\"hPa\",\n                             righttitlefontsize=16)\n\n# Set characteristics of text box\nprops = dict(facecolor='white', edgecolor='black', alpha=0.5)\n\n# Place text box\nax.text(0.40,\n        -0.1,\n        'CONTOUR FROM 948 TO 1064 BY 4',\n        transform=ax.transAxes,\n        fontsize=16,\n        bbox=props)\n\n# Make layout tight\nplt.tight_layout()\n\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
}