NBM_Probabilistic_Explorer.ipynb

{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "zVIP3cAxtGei"
      },
      "source": [
        "#**NBM Probabilistic Plotter**\n",
        "<a href=\"https://githubtocolab.com/csteele2/Wx4Colab/blob/master/NBM_Probabilistic_Explorer.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\"/></a> <br/>\n",
        "This notebook combines two preview notebooks to plot probailistic NBM temperature and precipitation data along with limited qualitative \"verification\" data (obs/LSRs). This notebook is meant to be used in the google colab environment and built with NWS SOOs in mind; so it may need adjustments if used elsewhere by others. Uses some code borrowed from [Brian Blaylock](https://github.com/blaylockbk) to more efficiently access the [NBM grib archive on AWS](https://noaa-nbm-grib2-pds.s3.amazonaws.com/index.html) for data, matplotlib, cartopy, and contextily to plot, [synoptic data preciptation api](https://developers.synopticdata.com/mesonet/v2/stations/precipitation/) for precipitation obs / [synoptic data](https://synopticdata.com/mesonet-api) statistics api for Min/Max T obs *(only tested with NOAA/NWS accounts)*, and [Iowa State LSR api](https://mesonet.agron.iastate.edu/request/gis/lsrs.phtml) for LSRs. Feel free to improve, steal, or use as a springboard for your own endeavours! <br/>\n",
        "<br/>\n",
        "Caleb Steele - https://github.com/csteele2/Wx4Colab\n",
        "<br/>\n",
        "\n",
        "> 21-Dec-2023: Just a few tweaks to keep it working with past events, and to also remove stamen maps.\n",
        "\n",
        "> 10-Oct-2023: Since cfgrib kinda stopped working with the existing implementation of condacolab, and since condacolab isn't necessary anymore now that google colab is on py 3.10+, removed the install of conda colab (which removed a step and error message! yay!)\n",
        "\n",
        "> 12-Sept-2023: Just a few tweaks and bug fixes. Changed the projection of CONUS from mercator so it isn't so distorted anymore (especially at higher latitudes). Added a workaround for the lowest (0.10\") threshold of psnow so it should also work again.\n",
        "\n",
        "> 12-Apr-2023: Added prob daily max wind plots. Since the precomputed threholds are more valid for old school marine forecasting, decided to add a function to interpolate to any desired threshold from the percentiles. Also applied this to PQPF, since I was lazy and didn't want to look up what thresholds are valid for 24, 48, and 72-hr PQPF (they are different). Unfortunately, since I switched to xarray/cfgrib, this needs a loop to decode the 99 percentiles individually, so it does add about 30 seconds to run time.\n",
        "\n",
        "> 17-Jan-2023: Some tweaks (like the CO domain) and updates to account for NBM v4.1 implementation. In NBM 4.1, PWPF (pSnow, pIce) are tied to PQPF, thus are not run every hour. However, they are not inserted into the QMD files, rather included in the next available core file, which is the QMD cycle + 7 hours (i.e. the 12Z QMD run PWPF is available in the 19Z core files, but PQPF, which it is tied to, retains the QMD cycling scheme. Which means the PQPF and PWPF from the same cycle will have different forecast hours). I account for this in the code, but its clear as mud, right?\n",
        "\n",
        "> 28-Dec-2022: A few bug fixes, especially related to PMxMnT option regarding legend, title, and filename\n",
        "\n",
        "> 22-Dec-2022: Initial release of combined notebook that uses xarray vice pygrib to read grib files"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "27yZ0aNmu1Up"
      },
      "source": [
        "## **1 - Install and Import**\n",
        "\n",
        "This will install and import everything we need. You need only run this once per session, then you can make all the changes to the form and make as many plots as you wish without having to rerun this cell (unless your session expires and spins down)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ChLVAGluX37I",
        "outputId": "cd36680c-2717-4a94-8e02-ee43eb5bfe8e"
      },
      "outputs": [],
      "source": [
        "!pip install ecmwflibs\n",
        "!pip install eccodes\n",
        "!pip install cfgrib\n",
        "!pip install cartopy contextily pyproj pyepsg netCDF4 xarray\n",
        "\n",
        "import numpy as np # to deal with arrays and math\n",
        "import xarray as xr # for reading gribs\n",
        "from urllib.request import urlretrieve # to get files\n",
        "import requests # to read a remote file\n",
        "import os, re, traceback # to do sorcery with subsetting nbm gribs\n",
        "\n",
        "import matplotlib # for plotting\n",
        "import matplotlib.patches as mpatches # helps us make a custom legend\n",
        "import matplotlib.pyplot as plt # to make actual plots\n",
        "import matplotlib.axes as maxes\n",
        "import matplotlib.patheffects as PathEffects # to add outline to text, etc.\n",
        "from matplotlib.path import Path\n",
        "from matplotlib.textpath import TextToPath\n",
        "from matplotlib.font_manager import FontProperties # we use this to use a custom font as markers\n",
        "matplotlib.rcParams['font.sans-serif'] = 'Liberation Sans'\n",
        "matplotlib.rcParams['font.family'] = \"sans-serif\"\n",
        "\n",
        "import pandas as pd # convenient to organize LSRs and obs\n",
        "import json # for parsing the syntopic api return\n",
        "from pyproj.crs import CRS # to warp maptiles\n",
        "\n",
        "from datetime import datetime, timedelta # to deal with datetime objects\n",
        "\n",
        "from cartopy import crs as ccrs, feature as cfeature # to make the plots maps\n",
        "from cartopy.io.shapereader import Reader # to read shapefiles\n",
        "from cartopy.feature import ShapelyFeature # also to read shapefiles\n",
        "import contextily as cx # to get fancy maptiles\n",
        "\n",
        "import warnings # to squash warnings that I have deemed insignificant\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3d1A0b4au8jV"
      },
      "source": [
        "## **2 - Edit Form Options & Go!**\n",
        "\n",
        "The cell below has some config things to set. For the curious, you can unhide the code hiding underneath."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 991
        },
        "id": "ykvnDFJNYUIB",
        "outputId": "555f3b39-c71c-4943-b021-96b400b51ff7"
      },
      "outputs": [],
      "source": [
        "from numpy.ma.core import flatten_mask\n",
        "#@title Plot Config { vertical-output: true, display-mode: \"form\" }\n",
        "# If this looks weird (like code instead of a form), try opening on Colab for the best experience!\n",
        "#@markdown ##<b>--------------- NBM Selection --------------- </b><br/>\n",
        "nbm_var = \"PQPF\" #@param [\"PQPF\", \"PMxMnT\", \"PMxWind\"]\n",
        "#@markdown Choose Valid (ending) date & hour (UTC) (hour not needed for PMxMnT or PMxWind)\n",
        "valid_date = \"2023-12-26\" #@param {type:\"date\"}\n",
        "valid_hour = 0 #@param {type:\"slider\", min:0, max:18, step:6}\n",
        "#@markdown Choose NBM QMD initialization date/hour (UTC) (PWPF is every hour, but PQPF is only 00/06/12/18Z).\n",
        "init_date = \"2023-12-20\" #@param {type:\"date\"}\n",
        "init_hour = 6 #@param {type:\"slider\", min:0, max:18, step:6}\n",
        "\n",
        "#@markdown ### <b>--------------- PQPF Options ---------------</b> <br/>\n",
        "#@markdown Choose the length of your forecast period (PQPF is only 24-hr pre-v4.1 upgrade (17-Jan-2023), otherwise 24, 48, or 72-hr available).\n",
        "forecast_length = \"48\" #@param [24, 48, 72]\n",
        "forecast_length = str(forecast_length)\n",
        "#@markdown Which precipitation type? (Snow, Ice, or QPF)\n",
        "ptype = \"snow\" #@param [\"snow\", \"ice\", \"qpf\"]\n",
        "#@markdown Do you want a probability of exceedance or a percentile?\n",
        "pqpf_var = \"prob\" #@param [\"prob\", \"percentile\"]\n",
        "#vt = \"2022-02-04T12:00\" #@param {type:\"string\"}\n",
        "#init = \"2022-01-31T13:00\" #@param {type:\"string\"}\n",
        "#@markdown Choose value or exceedance or percentile (you only need to change one / whichever variable you picked earlier)\n",
        "#@markdown Note: not all thresholds are available for all forecast lengths!\n",
        "threshold_qpf = 2 #@param type: number\n",
        "#threshold_qpf = \"{:.2f}\".format(threshold_qpf)\n",
        "#threshold_qpf = str(threshold_qpf)\n",
        "PQPF_thresh_dict = {\"24\":[\"0.01\", \"0.10\", \"0.25\", \"0.50\", \"1.00\", \"2.00\", \"3.00\", \"4.00\", \"5.00\", \"6.00\"],\n",
        "                    \"48\":[\"0.10\", \"1.00\", \"2.00\", \"4.00\", \"6.00\", \"8.00\", \"12.00\", \"18.00\", \"24.00\", \"30.00\"],\n",
        "                    \"72\":[\"0.10\", \"1.00\", \"2.00\", \"4.00\", \"6.00\", \"8.00\", \"12.00\", \"18.00\", \"24.00\", \"30.00\"]}\n",
        "#@markdown Do you want a tabbed overview? Otherwise, a single plot will be returned. (In development - only for pqpf for now)\n",
        "overview_opt = False#@param {type:\"boolean\"}\n",
        "\n",
        "#@markdown Otherwise, select threshold below.\n",
        "threshold_ice = \"0.10\" #@param [\"0.01\", \"0.10\", \"0.25\", \"0.50\", \"1.00\"]\n",
        "threshold_snow = \"1.00\" #@param [\"0.10\", \"1.00\", \"2.00\", \"4.00\", \"6.00\", \"8.00\", \"12.00\", \"18.00\", \"24.00\", \"30.00\"]\n",
        "pqpf_threshold_percentile = \"90\" #@param [\"5\", \"10\", \"25\", \"50\", \"75\", \"90\", \"95\"]\n",
        "threshold_dict = {\"0.01\":\"0.254\",\"0.10\":\"2.54\",\"0.25\":\"6.35\",\"0.50\":\"12.7\",\n",
        "                  \"1.00\":\"25.4\",\"2.00\":\"50.8\",\"3.00\":\"76.2\",\"4.00\":\"101.6\",\n",
        "                  \"5.00\":\"127\",\"6.00\":\"152.4\",\"8.00\":\"203.2\",\"12.00\":\"304.8\",\n",
        "                  \"18.00\":\"457.2\",\"24.00\":\"609.6\",\"30.00\":\"762\"}\n",
        "perc_suffix_dict = dict.fromkeys([\"5\",\"10\",\"25\",\"50\",\"75\",\"90\",\"95\", \"99\"],\"th\")\n",
        "perc_suffix_dict.update(dict.fromkeys(['1'], \"st\"))\n",
        "#if ptype == \"qpf\":\n",
        "  #threshold_mm = threshold_dict[threshold_qpf]\n",
        "if ptype == \"snow\":\n",
        "  threshold_mm = threshold_dict[threshold_snow]\n",
        "elif ptype == \"ice\":\n",
        "  threshold_mm = threshold_dict[threshold_ice]\n",
        "#threshold = \"025\" #@param {type:\"string\"}\n",
        "#@markdown <br />\n",
        "\n",
        "#@markdown ### <b>--------------- PMxMnT Options --------------- </b><br/>\n",
        "#@markdown Which temperature type? (MaxT or MinT)\n",
        "temp = \"MinT\" #@param [\"MinT\", \"MaxT\"]\n",
        "#@markdown Choose non/exeedance threshold\n",
        "#@markdown HotMaxT: 80, 90, 100, 110, 120 | ColdMaxT: 0, 28, 32 | ColdMinT: -40, -20,  0, 28, 32 | HotMinT: 80\n",
        "threshold_temp = 0 #@param [-40, -20, 0, 28, 32, 80, 90, 100, 110, 120]{type:\"raw\"}\n",
        "\n",
        "threshold_dict_FK = {-40:233,-20:244,0:255,28:270,\n",
        "                  32:273,80:299, 90:305,100:310,110:316,120:322}\n",
        "threshold_dict_prob = {\"5\":\"05\",\"10\":\"10\",\"25\":\"25\",\"50\":\"50\",\n",
        "                       \"75\":\"75\",\"90\":\"90\",\"95\":\"95\"}\n",
        "if temp == \"MaxT\":\n",
        "  if threshold_temp > 32:\n",
        "    element = \"HotMaxT\"\n",
        "  elif threshold_temp < 80:\n",
        "    element = \"ColdMaxT\"\n",
        "\n",
        "if temp == \"MinT\":\n",
        "  if threshold_temp > 32:\n",
        "    element = \"HotMinT\"\n",
        "  elif threshold_temp < 80:\n",
        "    element = \"ColdMinT\"\n",
        "#@markdown <br />\n",
        "\n",
        "#@markdown ### <b>--------------- PMxWind Options --------------- </b><br/>\n",
        "#@markdown Which type of wind: speed or gust?\n",
        "wind_type = \"gust\" #@param [\"speed\", \"gust\"]\n",
        "\n",
        "#@markdown Do you want a probability of exceedance or a percentile?\n",
        "wind_var = \"percentile\" #@param [\"prob\", \"percentile\"]\n",
        "\n",
        "#@markdown Choose value or exceedance or percentile (you only need to change one / whichever variable you picked earlier)\n",
        "wind_threshold = 39 #@param type: number\n",
        "wind_threshold_percentile = \"50\" #@param [\"5\", \"10\", \"25\", \"50\", \"75\", \"90\", \"95\"]\n",
        "\n",
        "\n",
        "#@markdown <br />\n",
        "\n",
        "#@markdown ##<b>--------------- Verification Selection --------------- </b><br/>\n",
        "#@markdown Do you want to plot LSRs?\n",
        "lsr_opt = False #@param {type:\"boolean\"}\n",
        "#@markdown Do you want to plot Obs? If you check this box, you also need to provide a synoptic token.\n",
        "ob_opt = False #@param {type:\"boolean\"}\n",
        "#@markdown Do you want to plot os misses?\n",
        "plot_ob_misses = False #@param {type:\"boolean\"}\n",
        "\n",
        "if nbm_var == \"PQPF\" and ptype != \"qpf\" and ob_opt:\n",
        "  print(\" ! > No obs option for anything other than QPF for PQPF. Fixing this for you, but next time, try to be more careful!\")\n",
        "  ob_opt = False\n",
        "#@markdown If obs selected, must paste in synoptic api token below - sign in or sign up at https://developers.synopticdata.com/signup/\n",
        "synoptic_token = \"\" #@param {type:\"string\"}\n",
        "#@markdown Which obs?\n",
        "network_selection = \"NWS\" #@param [\"NWS\", \"RAWS\", \"NWS+RAWS\", \"NWS+RAWS+HADS\", \"ALL\", \"CUSTOM\", \"LIST\"]\n",
        "#@markdown If Custom or List selected for network, enter comma separated network IDs (custom) or siteids (list)  WITH NO SPACES here. For help - https://developers.synopticdata.com/about/station-providers/\n",
        "network_input = \"106\"#@param {type:\"string\"}\n",
        "# Setup a dictionary for translating a form selection into a something we can pass to mesowest API\n",
        "network_dict = {\"NWS+RAWS+HADS\":\"&network=1,2,106\",\"NWS+RAWS\":\"&network=1,2\", \"NWS\":\"&network=1\", \"RAWS\": \"&network=2\", \"ALL\":\"\", \"CUSTOM\": \"&network=\"+network_input, \"LIST\": \"&stid=\"+network_input}\n",
        "network_string = network_dict[network_selection]\n",
        "\n",
        "#@markdown <br />\n",
        "\n",
        "#@markdown ##<b>--------------- Map Selection --------------- </b><br/>\n",
        "#@markdown Pick a map theme\n",
        "map_theme = \"Dark Grey Matter\" #@param [\"Light Shaded Relief\", \"Positron\", \"Dark Matter\", \"Dark Grey Matter\", \"ESRI Light Grey\", \"ESRI Dark Grey\"]\n",
        "if \"Dark\" in map_theme:\n",
        "  font_color='w'\n",
        "  face_color='#272727'\n",
        "else:\n",
        "  font_color = 'k'\n",
        "  face_color = 'w'\n",
        "#@markdown Set the map scale offset from default (i.e. a 1 would scale up one level [make labels bigger])\n",
        "map_scale_offset = \"-1\" #@param [\"-2\",\"-1\", \"0\", \"1\",\"2\"]\n",
        "map_zoom_offset = int(map_scale_offset)\n",
        "#@markdown Do you want CWA boundaries?\n",
        "cwa_opt = False #@param {type:\"boolean\"}\n",
        "#@markdown Pick your domain or select custom and input custom lat/lon\n",
        "dom = \"CONUS\" #@param [\"Custom\", \"CONUS\", \"ECONUS\", \"---- WESTERN REGION ----\", \"WR\",\"NR\",\"UT\",\"AZ\",\"SWUS\",\"PNW\", \"---- CENTRAL REGION ----\",\"CR\", \"NP\", \"GL\", \"CUS\", \"CO\", \"---- SOUTHERN REGION ----\", \"SR\", \"TXOK\", \"SE\", \"---- EASTERN REGION ----\", \"ER\", \"NE\"]\n",
        "#@markdown Enter custom lat/lon bounding box if custom was selected\n",
        "custom_name = \"NFR\" #@param {type:\"string\"}\n",
        "custom_bottom_lat =  39.8756#@param {type:\"number\"}\n",
        "custom_left_lon = -106.0893 #@param {type:\"number\"}\n",
        "custom_top_lat =  40.7665#@param {type:\"number\"}\n",
        "custom_right_lon = -104.5655 #@param {type:\"number\"}\n",
        "\n",
        "print(\"Figuring out what to do based on your inputs...\")\n",
        "\n",
        "derive_opt = False\n",
        "nbm_init = datetime.strptime(init_date,'%Y-%m-%d') + timedelta(hours=int(init_hour))\n",
        "nbm_init_string = nbm_init.strftime('%Y%m%d%H')\n",
        "if nbm_init < datetime(2023,1,17) and ptype == \"qpf\" and forecast_length != \"24\":\n",
        "  print(\" ❌ Ya goob! PQPF in v4.0 is only in 24-hr forecast length. I'll fix it for you this time, but next time try to me more careful!\")\n",
        "  forecast_length = \"24\"\n",
        "\n",
        "if nbm_var == \"PQPF\":\n",
        "  var = pqpf_var\n",
        "  threshold_percentile = pqpf_threshold_percentile\n",
        "  if var == \"prob\" and ptype == \"qpf\":\n",
        "    derive_opt = True\n",
        "elif nbm_var == \"PMxWind\":\n",
        "  var = wind_var\n",
        "  threshold_percentile = wind_threshold_percentile\n",
        "  forecast_length = 24\n",
        "  if wind_var == \"prob\":\n",
        "    derive_opt = True\n",
        "  ob_stat = \"maximum\"\n",
        "\n",
        "if nbm_var == \"PMxMnT\":\n",
        "  var = \"prob\"\n",
        "  if temp ==\"MaxT\":\n",
        "      nbm_qmd_valid_hour=\"06\"\n",
        "      valid_date = datetime.strptime(valid_date,'%Y-%m-%d')\n",
        "      valid_date_start = valid_date\n",
        "      valid_date_end = valid_date + timedelta(days=1)\n",
        "      obs_start_hour = \"1200\"\n",
        "      obs_end_hour = \"0600\"\n",
        "      ob_stat = \"maximum\"\n",
        "      valid_end_datetime = valid_date_end + timedelta(hours=(int(obs_end_hour)/100))\n",
        "      nbm_qmd_valid_end_datetime = valid_date_end + timedelta(hours=int(nbm_qmd_valid_hour))\n",
        "\n",
        "  elif temp == \"MinT\":\n",
        "      nbm_qmd_valid_hour=\"18\"\n",
        "      valid_date = datetime.strptime(valid_date,'%Y-%m-%d')\n",
        "      valid_date_start = valid_date\n",
        "      valid_date_end = valid_date\n",
        "      obs_start_hour = \"0000\"\n",
        "      obs_end_hour = \"1800\"\n",
        "      ob_stat = \"minimum\"\n",
        "      valid_end_datetime = valid_date_end + timedelta(hours=(int(obs_end_hour)/100))\n",
        "      nbm_qmd_valid_end_datetime = valid_date_end + timedelta(hours=int(nbm_qmd_valid_hour))\n",
        "\n",
        "elif nbm_var == \"PMxWind\":\n",
        "  nbm_qmd_valid_hour=\"06\"\n",
        "  valid_date = datetime.strptime(valid_date,'%Y-%m-%d')\n",
        "  valid_date_start = valid_date\n",
        "  valid_date_end = valid_date + timedelta(days=1)\n",
        "  obs_start_hour = \"0601\"\n",
        "  obs_end_hour = \"0600\"\n",
        "  ob_stat = \"maximum\"\n",
        "  valid_end_datetime = valid_date_end + timedelta(hours=(int(obs_end_hour)/100))\n",
        "  nbm_qmd_valid_end_datetime = valid_date_end + timedelta(hours=int(nbm_qmd_valid_hour))\n",
        "\n",
        "else:\n",
        "  nbm_qmd_valid_hour=str(valid_hour)\n",
        "  valid_date = datetime.strptime(valid_date,'%Y-%m-%d') + timedelta(hours=int(valid_hour))\n",
        "  # Below, we add a check to see if the NBM init date was after the implementation of 4.1.\n",
        "  # If so, PWPF and PQPF are tied together, but only PQPF is included in the QMD files,\n",
        "  # while PWPF is still left in the core files. However, there is a naming scheme difference\n",
        "  # we need to account for. QMD cycles are named after model init, NBM core cycles named after\n",
        "  # NBM run time (00Z QMD has 00Z models. 00Z NBM core is run at 00Z, and does not have 00Z models).\n",
        "  # So, add 7 hours to the NBM init time to account for this (PWPF only)\n",
        "  if ptype != \"qpf\":\n",
        "    #if nbm_init >= datetime(2023,1,17):\n",
        "    nbm_init = nbm_init + timedelta(hours=7)\n",
        "    nbm_init_string = nbm_init.strftime('%Y%m%d%H')\n",
        "  valid_date_start = valid_date - timedelta(hours=int(forecast_length))\n",
        "  valid_date_end = valid_date\n",
        "  obs_start_hour = str(valid_hour)+\"00\"\n",
        "  obs_end_hour = str(valid_hour)+\"00\"\n",
        "  ob_stat = \"total\"\n",
        "  nbm_qmd_valid_end_datetime = valid_date_end\n",
        "\n",
        "nbm_qmd_fhdelta = nbm_qmd_valid_end_datetime - nbm_init\n",
        "nbm_qmd_forecasthour = nbm_qmd_fhdelta.total_seconds() / 3600.\n",
        "if nbm_var == \"PMxMnT\":\n",
        "  nbm_qmd_forecasthour_start = nbm_qmd_forecasthour - 18\n",
        "else:\n",
        "  nbm_qmd_forecasthour_start = nbm_qmd_forecasthour - int(forecast_length)\n",
        "\n",
        "if nbm_qmd_forecasthour < 0:\n",
        "  print(\" ❌ Uh oh! Your NBM initialization is AFTER your desired valid date. Did you get those backwards?\")\n",
        "\n",
        "\n",
        "def timestring(modeltime):\n",
        "    outtime = '%s' % (modeltime.strftime('%HZ %a %m-%d-%Y'))\n",
        "    return outtime\n",
        "\n",
        "def timestring2(modeltime):\n",
        "  outtime = '%s' % (modeltime.strftime('%HZ %m-%d-%Y'))\n",
        "  return outtime\n",
        "\n",
        "initime = timestring2(nbm_init)\n",
        "CURTIMESTRING = timestring(valid_date)\n",
        "\n",
        "latloncrs = ccrs.PlateCarree()\n",
        "#proj = ccrs.epsg(3857)\n",
        "if dom == \"CONUS\":\n",
        "  proj = ccrs.AlbersEqualArea(central_longitude=-98.35,\n",
        "                              central_latitude=39.5,\n",
        "                              false_easting=0.0,\n",
        "                              false_northing=0.0,\n",
        "                              standard_parallels=(20.0, 50.0),\n",
        "                              globe=None)\n",
        "else:\n",
        "  proj = ccrs.Mercator.GOOGLE\n",
        "\n",
        "width = 7 # sets figure width value\n",
        "height = 7 # sets figure height value\n",
        "\n",
        "domain_dict = {\"Custom\":{\"west\":custom_left_lon,\n",
        "                        \"east\":custom_right_lon,\n",
        "                        \"north\":custom_top_lat,\n",
        "                        \"south\":custom_bottom_lat,\n",
        "                        \"zoom_adj\": 1,\n",
        "                        \"legend\":3},\n",
        "\n",
        "               \"CONUS\":{\"west\":-123.650,\n",
        "                    \"south\":23.377,\n",
        "                    \"east\":-71.488,\n",
        "                    \"north\":50.924,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"ECONUS\":{\"west\":-104.36,\n",
        "                    \"south\":24.735,\n",
        "                    \"east\":-66.453,\n",
        "                    \"north\":49.755,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"WR\":{\"west\":-126.917,\n",
        "                    \"south\":30.586,\n",
        "                    \"east\":-102.740,\n",
        "                    \"north\":49.755,\n",
        "                    \"zoom_adj\": 1,\n",
        "                     \"legend\":4},\n",
        "\n",
        "               \"UT\":{\"west\":-117.02,\n",
        "                      \"east\":-106.92,\n",
        "                      \"north\":42.13,\n",
        "                      \"south\":36.80,\n",
        "                      \"zoom_adj\": 1,\n",
        "                     \"legend\":4},\n",
        "\n",
        "               \"NR\":{\"west\":-117.5177,\n",
        "                    \"south\":41.9071,\n",
        "                    \"east\":-103.38071,\n",
        "                    \"north\":49.3085,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"PNW\":{\"west\":-125.4510,\n",
        "                    \"south\":41.8754,\n",
        "                    \"east\":-110.9318,\n",
        "                    \"north\":49.5767,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"SWUS\":{\"west\":-125.582,\n",
        "                    \"south\":31.136,\n",
        "                    \"east\":-108.689,\n",
        "                    \"north\":42.859,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":3},\n",
        "\n",
        "               \"AZ\":{\"west\":-115.596,\n",
        "                    \"south\":31.113,\n",
        "                    \"east\":-107.887,\n",
        "                    \"north\":37.446,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":4},\n",
        "\n",
        "\n",
        "               \"CR\":{\"west\":-111.534,\n",
        "                    \"south\":35.118,\n",
        "                    \"east\":-82.263,\n",
        "                    \"north\":49.755,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"NP\":{\"west\":-105.244,\n",
        "                    \"south\":42.173,\n",
        "                    \"east\":-89.426,\n",
        "                    \"north\":49.474,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"GL\":{\"west\":-97.606,\n",
        "                    \"south\":38.735,\n",
        "                    \"east\":-74.916,\n",
        "                    \"north\":49.292,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":3},\n",
        "\n",
        "               \"CUS\":{\"west\":-111.553,\n",
        "                    \"south\":34.794,\n",
        "                    \"east\":-88.533,\n",
        "                    \"north\":46.357,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":3},\n",
        "\n",
        "               \"CO\":{\"west\":-110.42426,\n",
        "                    \"south\":35.93570,\n",
        "                    \"east\":-100.80419,\n",
        "                    \"north\":41.85529,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"SR\":{\"west\":-109.758,\n",
        "                    \"south\":23.313,\n",
        "                    \"east\":-78.247,\n",
        "                    \"north\":37.899,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":3},\n",
        "\n",
        "               \"TXOK\":{\"west\":-106.95,\n",
        "                    \"south\":26.06,\n",
        "                    \"east\":-86.76,\n",
        "                    \"north\":37.76,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"SE\":{\"west\":-92.974,\n",
        "                    \"south\":24.578,\n",
        "                    \"east\":-75.1311,\n",
        "                    \"north\":37.390,\n",
        "                    \"zoom_adj\": 1,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"ER\":{\"west\":-85.629,\n",
        "                    \"south\":31.723,\n",
        "                    \"east\":-66.465,\n",
        "                    \"north\":47.676,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":4},\n",
        "\n",
        "               \"NE\":{\"west\":-85.629,\n",
        "                    \"south\":37.654,\n",
        "                    \"east\":-66.00,\n",
        "                    \"north\":47.825,\n",
        "                    \"zoom_adj\": 0,\n",
        "                    \"legend\":4},\n",
        "}\n",
        "\n",
        "west = domain_dict[dom][\"west\"]\n",
        "south = domain_dict[dom][\"south\"]\n",
        "east = domain_dict[dom][\"east\"]\n",
        "north = domain_dict[dom][\"north\"]\n",
        "#map_zoom_offset = domain_dict[dom][\"zoom_adj\"]\n",
        "LLOC = domain_dict[dom][\"legend\"]\n",
        "\n",
        "icon_dict={\"HotMaxT\":\n",
        "            {\"Ob_hit\":\n",
        "              {\"icon\":\"\\ue040\",\"color\":\"#4f0615\",\"label\":f'Ob > {threshold_temp}$^\\circ$F'},\n",
        "             \"Ob_miss\":\n",
        "              {\"icon\":\"\\ue03f\",\"color\":\"#bfbfbf\",\"label\":f'Ob ≤ {threshold_temp}$^\\circ$F'},\n",
        "            },\n",
        "           \"ColdMinT\":\n",
        "            {\"Ob_hit\":\n",
        "              {\"icon\":\"\\ue03f\",\"color\":\"navy\",\"label\":f'Ob < {threshold_temp}$^\\circ$F'},\n",
        "             \"Ob_miss\":\n",
        "              {\"icon\":\"\\ue040\",\"color\":\"#bfbfbf\",\"label\":f'Ob ≥ {threshold_temp}$^\\circ$F'},\n",
        "            },\n",
        "           \"ColdMaxT\":\n",
        "            {\"Ob_hit\":\n",
        "              {\"icon\":\"\\ue03f\",\"color\":\"navy\",\"label\":f'Ob < {threshold_temp}$^\\circ$F'},\n",
        "             \"Ob_miss\":\n",
        "              {\"icon\":\"\\ue040\",\"color\":\"#bfbfbf\",\"label\":f'Ob ≥ {threshold_temp}$^\\circ$F'},\n",
        "            },\n",
        "           \"HotMinT\":\n",
        "            {\"Ob_hit\":\n",
        "              {\"icon\":\"\\ue040\",\"color\":\"#4f0615\",\"label\":f'Ob > {threshold_temp}$^\\circ$F'},\n",
        "             \"Ob_miss\":\n",
        "              {\"icon\":\"\\ue03f\",\"color\":\"#bfbfbf\",\"label\":f'Ob ≤ {threshold_temp}$^\\circ$F'},\n",
        "            },\n",
        "           \"qpf\":\n",
        "            {\"LSR\":\n",
        "              {\"heavyrain\":{\"icon\":\"\\uf740\",\"color\":\"tab:orange\",\"label\":str(float(threshold_qpf))+\"\\\"+ LSR\"},\n",
        "              \"flooding\":{\"icon\":\"\\uf773\",\"color\":\"tab:red\",\"label\":\"Flood LSR\"},\n",
        "            },\n",
        "             \"Ob\":\n",
        "              {\"icon\":\"\\uf043\",\"color\":\"#033751\",\"label\":str(float(threshold_qpf))+\"\\\"+ Ob\"},\n",
        "            },\n",
        "           \"snow\":\n",
        "            {\"LSR\":\n",
        "              {\"icon\":\"\\uf2dc\",\"color\":\"#E98FAC\",\"label\":str(float(threshold_snow))+\"\\\"+ LSR\"},\n",
        "             },\n",
        "           \"ice\":\n",
        "            {\"LSR\":\n",
        "              {\"icon\":\"\\uf7ad\",\"color\":\"#EC8E1C\",\"label\":str(float(threshold_ice))+\"\\\"+ LSR\"},\n",
        "             },\n",
        "           \"PMxWind\":\n",
        "            {\"LSR\":\n",
        "              {\"icon\":\"\\uf72e\",\"color\":\"#FFC000\",\"label\":str(int(wind_threshold))+\"+ MPH LSR\"},\n",
        "             \"Ob_hit\":\n",
        "              {\"icon\":\"\\ue040\",\"color\":\"#FFC000\",\"label\":f'Ob > {wind_threshold} MPH'},\n",
        "             \"Ob_miss\":\n",
        "              {\"icon\":\"\\ue03f\",\"color\":\"#bfbfbf\",\"label\":f'Ob ≤ {wind_threshold} MPH'},\n",
        "             },\n",
        "}\n",
        "\n",
        "\n",
        "\n",
        "########################################################################################################################\n",
        "# Common functions here                                                                                                #\n",
        "########################################################################################################################\n",
        "\n",
        "# The following is a function to grap a font and use the glyphs as markers in a scatter plot\n",
        "if os.path.exists(\"Font_Awesome_6_Free-Solid-900.otf\"):\n",
        "  pass\n",
        "else:\n",
        "  urlretrieve(\"https://www.dynamicmeteorology.com/style/Font_Awesome_6_Free-Solid-900.otf\", \"Font_Awesome_6_Free-Solid-900.otf\")\n",
        "fp = FontProperties(fname=r\"Font_Awesome_6_Free-Solid-900.otf\")\n",
        "\n",
        "def get_marker(symbol):\n",
        "  v, codes = TextToPath().get_text_path(fp, symbol)\n",
        "  v = np.array(v)\n",
        "  mean = np.mean([np.max(v,axis=0), np.min(v, axis=0)], axis=0)\n",
        "  return Path(v-mean, codes, closed=False)\n",
        "\n",
        "def mm_to_in(millimeters):\n",
        "  inches = millimeters / 25.4\n",
        "  return inches\n",
        "\n",
        "def m_to_in(meters):\n",
        "  inches = meters * 39.3701\n",
        "  return inches\n",
        "\n",
        "def K_to_F(kelvin):\n",
        "  fahrenheit = 1.8*(kelvin-273)+32.\n",
        "  return fahrenheit\n",
        "\n",
        "def ms_to_mph(ms):\n",
        "  mph = ms * 2.237\n",
        "  return mph\n",
        "\n",
        "def find_roots(x,y):\n",
        "  s = np.abs(np.diff(np.sign(y))).astype(bool)\n",
        "  return x[:-1][s] + np.diff(x)[s]/(np.abs(y[1:][s]/y[:-1][s])+1)\n",
        "\n",
        "def compute_exceedance_probability(threshold, percentile_grids):\n",
        "  # Must have percentile grids as 1st-99th\n",
        "  # Find the index of the percentile grid containing the threshold value\n",
        "  percentile_grids2 = np.copy(percentile_grids)\n",
        "  percentile_grids[percentile_grids > threshold] = 0\n",
        "  threshold_percentile = np.argmax(percentile_grids, axis=0)\n",
        "  threshold_percentile = np.reshape(threshold_percentile, percentile_grids[0].shape)\n",
        "  probability = (100 - (threshold_percentile+1))\n",
        "  probability[percentile_grids2[98,:,:] < threshold] = 0\n",
        "\n",
        "  return probability\n",
        "\n",
        "\n",
        "# This bit of code to subset the grib and only download things we want is\n",
        "# shamelessly stolen from Brian Blaylock (https://github.com/blaylockbk)\n",
        "def download_subset(remote_url, remote_file, local_filename):\n",
        "  print(\"  > Downloading a subset of NBM gribs\")\n",
        "  local_file = \"nbm/\"+local_filename\n",
        "  if nbm_var == \"PQPF\":\n",
        "    if \"qmd\" in remote_file:\n",
        "      if ptype == \"qpf\":\n",
        "        if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "          search_string = f':APCP:surface:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day acc fcst:'\n",
        "        else:\n",
        "          search_string = f':APCP:surface:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour acc fcst:'\n",
        "    if \"core\" in remote_file:\n",
        "      if ptype == \"snow\":\n",
        "        if nbm_init > datetime(2022,7,5):\n",
        "          if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "            search_string = f':ASNOW:surface:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day acc'\n",
        "          else:\n",
        "            search_string = f':ASNOW:surface:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour acc'\n",
        "        else:\n",
        "          search_string = f':ASNOW:surface:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour acc'\n",
        "      elif ptype == \"ice\":\n",
        "        if nbm_init > datetime(2022,7,5):\n",
        "          if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "            search_string = f':FICEAC:surface:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day acc'\n",
        "          else:\n",
        "            search_string = f':FICEAC:surface:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour acc'\n",
        "        else:\n",
        "          search_string = f'228:surface:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour acc'\n",
        "  elif nbm_var == \"PMxMnT\":\n",
        "    if temp == \"MaxT\":\n",
        "      if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "        search_string = f':TMP:2 m above ground:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day max fcst:prob'\n",
        "      else:\n",
        "        search_string = f':TMP:2 m above ground:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour max fcst:prob'\n",
        "    elif temp == \"MinT\":\n",
        "      if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "        search_string = f':TMP:2 m above ground:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day min fcst:prob'\n",
        "      else:\n",
        "        search_string = f':TMP:2 m above ground:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour min fcst:prob'\n",
        "  elif nbm_var == \"PMxWind\":\n",
        "    if wind_type == \"speed\":\n",
        "      if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "        search_string = f':WIND:10 m above ground:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day max fcst:'\n",
        "      else:\n",
        "        search_string = f':WIND:10 m above ground:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour max fcst:'\n",
        "    elif wind_type == \"gust\":\n",
        "      if (int(nbm_qmd_forecasthour_start) % 24 == 0) and (int(nbm_qmd_forecasthour) % 24 ==0):\n",
        "        search_string = f':GUST:10 m above ground:{str(int(int(nbm_qmd_forecasthour_start)/24))}-{str(int(int(nbm_qmd_forecasthour)/24))} day max fcst:'\n",
        "      else:\n",
        "        search_string = f':GUST:10 m above ground:{str(int(nbm_qmd_forecasthour_start))}-{str(int(nbm_qmd_forecasthour))} hour max fcst:'\n",
        "\n",
        "  #print(search_string)\n",
        "  idx = remote_url+\".idx\"\n",
        "  r = requests.get(idx)\n",
        "  if not r.ok:\n",
        "    print('     ❌ SORRY! Status Code:', r.status_code, r.reason)\n",
        "    print(f'      ❌ It does not look like the index file exists: {idx}')\n",
        "\n",
        "  lines = r.text.split('\\n')\n",
        "  expr = re.compile(search_string)\n",
        "  byte_ranges = {}\n",
        "  for n, line in enumerate(lines, start=1):\n",
        "      # n is the line number (starting from 1) so that when we call for\n",
        "      # `lines[n]` it will give us the next line. (Clear as mud??)\n",
        "\n",
        "      # Use the compiled regular expression to search the line\n",
        "      if expr.search(line):\n",
        "          # aka, if the line contains the string we are looking for...\n",
        "\n",
        "          # Get the beginning byte in the line we found\n",
        "          parts = line.split(':')\n",
        "          rangestart = int(parts[1])\n",
        "\n",
        "          # Get the beginning byte in the next line...\n",
        "          if n+1 < len(lines):\n",
        "              # ...if there is a next line\n",
        "              parts = lines[n].split(':')\n",
        "              rangeend = int(parts[1])\n",
        "          else:\n",
        "              # ...if there isn't a next line, then go to the end of the file.\n",
        "              rangeend = ''\n",
        "\n",
        "          # Store the byte-range string in our dictionary,\n",
        "          # and keep the line information too so we can refer back to it.\n",
        "          byte_ranges[f'{rangestart}-{rangeend}'] = line\n",
        "          #print(line)\n",
        "  for i, (byteRange, line) in enumerate(byte_ranges.items()):\n",
        "\n",
        "        if i == 0:\n",
        "            # If we are working on the first item, overwrite the existing file.\n",
        "            curl = f'curl -s --range {byteRange} {remote_url} > {local_file}'\n",
        "        else:\n",
        "            # If we are working on not the first item, append the existing file.\n",
        "            curl = f'curl -s --range {byteRange} {remote_url} >> {local_file}'\n",
        "        try:\n",
        "          if nbm_var == \"PQPF\" or nbm_var == \"PMxWind\":\n",
        "            num, byte, date, var, level, forecast, _ = line.split(':')\n",
        "          elif nbm_var == \"PMxMnT\":\n",
        "            num, byte, date, var, level, forecast, fthresh, ftype = line.split(':')\n",
        "        except:\n",
        "          pass\n",
        "\n",
        "        #print(f'  Downloading GRIB line [{num:>3}]: variable={var}, level={level}, forecast={forecast}')\n",
        "        os.system(curl)\n",
        "\n",
        "  if os.path.exists(local_file):\n",
        "      print(f'      ✅ Success! Searched for [{search_string}] and got [{len(byte_ranges)}] GRIB fields and saved as {local_file}')\n",
        "      return local_file\n",
        "  else:\n",
        "      print(print(f'      ❌ Unsuccessful! Searched for [{search_string}] and did not find anything!'))\n",
        "\n",
        "\n",
        "def get_nbm():\n",
        "  print('Getting and processing NBM...')\n",
        "  global nbm, nbmlats, nbmlons\n",
        "\n",
        "  percList = np.arange(1,100,1)\n",
        "\n",
        "  if os.path.exists(\"nbm\"):\n",
        "    pass\n",
        "  else:\n",
        "    os.system('mkdir nbm')\n",
        "\n",
        "  nbm_init_filen = nbm_init.strftime('%Y%m%d') + \"_\" + nbm_init.strftime('%H')\n",
        "  nbm_url_base = \"https://noaa-nbm-grib2-pds.s3.amazonaws.com/blend.\"+nbm_init.strftime('%Y%m%d') \\\n",
        "              +\"/\"+nbm_init.strftime('%H')+\"/\"\n",
        "  if nbm_var == \"PQPF\" and ptype == \"qpf\":\n",
        "    nbm_file = f\"blend.t{int(nbm_init.strftime('%H')):02}z.qmd.f{int(nbm_qmd_forecasthour):03}.co.grib2\"\n",
        "    nbm_file_subset = f\"blend.t{int(nbm_init.strftime('%H')):02}z.qmd.{nbm_init_string}f{int(nbm_qmd_forecasthour):03}.co.pqpf{str(forecast_length)}_subset.grib2\"\n",
        "    nbm_url = nbm_url_base+\"qmd/\"+nbm_file\n",
        "  elif nbm_var == \"PQPF\" and ptype != \"qpf\":\n",
        "    nbm_file = f\"blend.t{int(nbm_init.strftime('%H')):02}z.core.f{int(nbm_qmd_forecasthour):03}.co.grib2\"\n",
        "    nbm_file_subset = f\"blend.t{int(nbm_init.strftime('%H')):02}z.core.{nbm_init_string}f{int(nbm_qmd_forecasthour):03}.co.p{ptype}{str(forecast_length)}_subset.grib2\"\n",
        "    nbm_url = nbm_url_base+\"core/\"+nbm_file\n",
        "  elif nbm_var == \"PMxMnT\":\n",
        "    nbm_file = f\"blend.t{int(nbm_init.strftime('%H')):02}z.qmd.f{int(nbm_qmd_forecasthour):03}.co.grib2\"\n",
        "    nbm_file_subset = f\"blend.t{int(nbm_init.strftime('%H')):02}z.qmd.{nbm_init_filen}f{int(nbm_qmd_forecasthour):03}.co.{element}_subset.grib2\"\n",
        "    nbm_url = nbm_url_base+\"qmd/\"+nbm_file\n",
        "  elif nbm_var == \"PMxWind\":\n",
        "    nbm_file = f\"blend.t{int(nbm_init.strftime('%H')):02}z.qmd.f{int(nbm_qmd_forecasthour):03}.co.grib2\"\n",
        "    nbm_file_subset = f\"blend.t{int(nbm_init.strftime('%H')):02}z.qmd.{nbm_init_string}f{int(nbm_qmd_forecasthour):03}.co.p{wind_type}_subset.grib2\"\n",
        "    nbm_url = nbm_url_base+\"qmd/\"+nbm_file\n",
        "\n",
        "  if os.path.exists(\"nbm/\"+nbm_file_subset):\n",
        "    print(\"  > NBM probabilistic already exists\")\n",
        "  else:\n",
        "    download_subset(nbm_url, nbm_file, nbm_file_subset)\n",
        "\n",
        "  if derive_opt == False:\n",
        "    if nbm_var == \"PMxMnT\":\n",
        "      if \"MaxT\" in element:\n",
        "        if \"Hot\" in element:\n",
        "          cf_args = {'filter_by_keys':{'probabilityTypeName':\"Probability of event above upper limit\",'upperLimit':threshold_dict_FK[threshold_temp]}}\n",
        "        elif \"Cold\" in element:\n",
        "          cf_args = {'filter_by_keys':{'probabilityTypeName':\"Probability of event below lower limit\",'lowerLimit':threshold_dict_FK[threshold_temp]}}\n",
        "      elif \"MinT\" in element:\n",
        "        if \"Cold\" in element:\n",
        "          cf_args = {'filter_by_keys':{'probabilityTypeName':\"Probability of event below lower limit\",'lowerLimit':threshold_dict_FK[threshold_temp]}}\n",
        "        elif \"Hot\" in element:\n",
        "          cf_args = {'filter_by_keys':{'probabilityTypeName':\"Probability of event above upper limit\",'upperLimit':threshold_dict_FK[threshold_temp]}}\n",
        "    elif nbm_var == \"PQPF\":\n",
        "      if var == \"percentile\":\n",
        "        cf_args = {'filter_by_keys':{'percentileValue':int(threshold_percentile)}}\n",
        "      elif var == \"prob\":\n",
        "        if ptype !=\"snow\":\n",
        "          cf_args = {'filter_by_keys':{'upperLimit':float(threshold_mm), 'probabilityTypeName':\"Probability of event above upper limit\"}}\n",
        "        else:\n",
        "          # Not sure why 0.1 inches is encoded as 0.0025399999999999997 meters instead of 0.00254 now\n",
        "          if ptype == \"snow\" and threshold_snow == \"0.10\":\n",
        "            thresholdKey = 0.0025399999999999997\n",
        "          else:\n",
        "            thresholdKey = float(threshold_mm)/1000.\n",
        "          cf_args = {'filter_by_keys':{'upperLimit':thresholdKey, 'probabilityTypeName':\"Probability of event above upper limit\"}}\n",
        "    elif nbm_var == \"PMxWind\":\n",
        "      cf_args = {'filter_by_keys':{'percentileValue':int(threshold_percentile)}}\n",
        "\n",
        "    with xr.open_dataset(\"nbm/\"+nbm_file_subset, engine='cfgrib', backend_kwargs=cf_args) as nbmgrb:\n",
        "      #print(nbmgrb)\n",
        "      nbmlats, nbmlons = nbmgrb.latitude.values, nbmgrb.longitude.values\n",
        "      datakey = list(nbmgrb.data_vars.keys())[0]\n",
        "      if nbm_var == \"PMxMnT\":\n",
        "        nbm = nbmgrb[datakey].values\n",
        "      elif nbm_var == \"PQPF\" and ptype == \"snow\":\n",
        "        if var == \"percentile\":\n",
        "          nbm = m_to_in(nbmgrb[datakey].values)\n",
        "        else:\n",
        "          nbm = nbmgrb[datakey].values\n",
        "      elif nbm_var == \"PQPF\" and ptype == \"ice\":\n",
        "        if var == \"percentile\":\n",
        "          nbm = mm_to_in(nbmgrb[datakey].values)\n",
        "        else:\n",
        "          nbm = nbmgrb[datakey].values\n",
        "      elif nbm_var == \"PQPF\" and ptype == \"qpf\":\n",
        "        if var == \"percentile\":\n",
        "          nbm = mm_to_in(nbmgrb[datakey].values)\n",
        "        else:\n",
        "          nbm = nbmgrb[datakey].values\n",
        "      elif nbm_var == \"PMxWind\" and var == \"percentile\":\n",
        "        if wind_type == \"speed\":\n",
        "          nbm = ms_to_mph(nbmgrb[datakey].values)\n",
        "        elif wind_type == \"gust\":\n",
        "          nbm = ms_to_mph(nbmgrb[datakey].values)\n",
        "\n",
        "  elif derive_opt == True:\n",
        "    print(\"  > Deriving probability of exceedance\")\n",
        "    percentile_cube = []\n",
        "    if nbm_var == \"PMxWind\":\n",
        "      threshold = wind_threshold\n",
        "      for percentile in percList:\n",
        "        with xr.open_dataset(\"nbm/\"+nbm_file_subset, engine='cfgrib', backend_kwargs={'filter_by_keys':{'percentileValue':int(percentile)}}) as nbmgrb:\n",
        "          datakey = list(nbmgrb.data_vars.keys())[0]\n",
        "          if percentile == 1:\n",
        "            nbmlats, nbmlons = nbmgrb.latitude.values, nbmgrb.longitude.values\n",
        "          if wind_type == \"gust\":\n",
        "            percentile_cube.append(ms_to_mph(nbmgrb[datakey].values))\n",
        "          elif wind_type == \"speed\":\n",
        "            percentile_cube.append(ms_to_mph(nbmgrb[datakey].values))\n",
        "    elif nbm_var == \"PQPF\" and ptype == \"qpf\":\n",
        "      threshold = threshold_qpf\n",
        "      for percentile in percList:\n",
        "        with xr.open_dataset(\"nbm/\"+nbm_file_subset, engine='cfgrib', backend_kwargs={'filter_by_keys':{'percentileValue':int(percentile)}}) as nbmgrb:\n",
        "          datakey = list(nbmgrb.data_vars.keys())[0]\n",
        "          if percentile == 1:\n",
        "            nbmlats, nbmlons = nbmgrb.latitude.values, nbmgrb.longitude.values\n",
        "          percentile_cube.append(mm_to_in(nbmgrb[datakey].values))\n",
        "\n",
        "    #threshold_grid = np.ones_like(percentile_cube[0]) * threshold\n",
        "    percentile_cube = np.array(percentile_cube)\n",
        "    nbm = compute_exceedance_probability(threshold, percentile_cube)\n",
        "\n",
        "    print(f'   ✅ Got it! Max: {int(np.max(nbm))}')\n",
        "\n",
        "\n",
        "\n",
        "### set up the blank plot with mapstuffs\n",
        "def blankmap():\n",
        "  global maplayertext1, maplayertext2\n",
        "  print(\" > Initializing map\")\n",
        "  plt.figure(figsize=(width,height),frameon=True, facecolor=face_color)\n",
        "  F = plt.gcf()  # Gets the current figure\n",
        "  ax = plt.axes(projection=proj)\n",
        "\n",
        "### Here is where you set up the domain.\n",
        "### Want to add another? Just copy the last (elif) one and change the bounds (try to keep it square)\n",
        "### Note the attributes are turned OFF on the cx.add_basemap layers IF you have mixed and matched provider sources. \\\n",
        "### This is because each attribution goes on top of the other, and thus are manually added so they remain legible.\n",
        "  zoom = (cx.tile._calculate_zoom(west, south, east, north) - map_zoom_offset)\n",
        "\n",
        "  ax.set_extent([west, east, south, north], crs=latloncrs)\n",
        "\n",
        "  print(' > Adding fancy map tiles')\n",
        "  maplayertext2 = \"\"\n",
        "  maylayertext1 = \"\"\n",
        "  if map_theme == \"Light Shaded Relief\":\n",
        "    ax.add_feature(cfeature.LAND, edgecolor='none', facecolor='#FAFAF8', zorder=-2)\n",
        "    cx.add_basemap(ax, source=\"https://basemap.nationalmap.gov/arcgis/rest/services/USGSShadedReliefOnly/MapServer/tile/{z}/{y}/{x}\",\n",
        "                    attribution=False, crs=proj, alpha =0.8, zorder=-1)\n",
        "    maplayertext1 = \"Map tiles: © USGS Earth Resources Observation & Science (EROS) Center: GMTED2010\" #for USGS Hillshade\n",
        "    ax.add_feature(cfeature.OCEAN, edgecolor='none', facecolor='#b3bbbd', zorder=1) # adds fill over the ocean\n",
        "    ax.add_feature(cfeature.LAKES, edgecolor='none', facecolor='#b3bbbd', zorder=1) # adds fill over lakes\n",
        "    cx.add_basemap(ax, source=\"http://services.arcgisonline.com/arcgis/rest/services/Reference/World_Boundaries_and_Places/MapServer/tile/{z}/{y}/{x}\",\n",
        "                    attribution=False, crs=proj, zoom=zoom, zorder=50)\n",
        "    cx.add_basemap(ax, source='http://services.arcgisonline.com/arcgis/rest/services/Reference/World_Transportation/MapServer/tile/{z}/{y}/{x}',\n",
        "                    crs=proj, zoom=zoom, zorder=49)\n",
        "    maplayertext2 = \"Esri, HERE, Garmin, OpenStreetMap contributors\"\n",
        "\n",
        "  elif map_theme == \"Positron\":\n",
        "    cx.add_basemap(ax, source=cx.providers.CartoDB.PositronNoLabels, zorder=-1, attribution=False, crs=proj)\n",
        "    ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', edgecolor='#e5e2e3', facecolor='none', linewidth=0.5, zorder=48, alpha=0.5)) # adds state borders\n",
        "    ax.add_feature(cfeature.COASTLINE, edgecolor=\"#FAFAF8\", facecolor=None, linewidth=0.25, alpha=0.5, zorder=48) # adds coastline (since that isn't included)\n",
        "    ax.add_feature(cfeature.BORDERS, edgecolor=\"#FAFAF8\", facecolor=None, linewidth=0.5, alpha=0.75, zorder=48) # adds country borders (since that isn't included)\n",
        "    cx.add_basemap(ax, source=cx.providers.CartoDB.PositronOnlyLabels, zoom=zoom, zorder=49, attribution=False, crs=proj)\n",
        "    maplayertext1 = \"Map tiles: \" + str(cx.providers.CartoDB.Positron.attribution)\n",
        "\n",
        "  elif map_theme == \"Dark Grey Matter\":\n",
        "    ax.add_feature(cfeature.LAND, edgecolor='none', facecolor='#414143', zorder=-2)\n",
        "    ax.add_feature(cfeature.OCEAN, edgecolor='none', facecolor='#232227', zorder=1) # adds fill over the ocean\n",
        "    ax.add_feature(cfeature.LAKES, edgecolor='none', facecolor='#232227', zorder=1) # adds fill over lakes\n",
        "    #cx.add_basemap(ax, source=cx.providers.CartoDB.DarkMatterNoLabels, zorder=-1, attribution=False)\n",
        "    ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', edgecolor='#1c1c1c', facecolor='none', linewidth=0.5, zorder=48, alpha=0.5)) # adds state borders\n",
        "    ax.add_feature(cfeature.COASTLINE, edgecolor='k', facecolor=None, linewidth=0.25, alpha=0.5, zorder=48) # adds coastline (since that isn't included)\n",
        "    ax.add_feature(cfeature.BORDERS, edgecolor='k', facecolor=None, linewidth=0.5, alpha=0.75, zorder=48) # adds country borders (since that isn't included)\n",
        "    cx.add_basemap(ax, source=cx.providers.CartoDB.DarkMatterOnlyLabels, zoom=zoom, zorder=49, attribution=False, crs=proj)\n",
        "    maplayertext1 = \"Map tiles: \" + str(cx.providers.CartoDB.DarkMatter.attribution)\n",
        "\n",
        "  elif map_theme == \"Dark Matter\":\n",
        "    cx.add_basemap(ax, source=cx.providers.CartoDB.DarkMatterNoLabels, zorder=-1, attribution=False, crs=proj)\n",
        "    ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', edgecolor='#1c1c1c', facecolor='none', linewidth=0.5, zorder=48, alpha=0.5)) # adds state borders\n",
        "    ax.add_feature(cfeature.COASTLINE, edgecolor='k', facecolor=None, linewidth=0.25, alpha=0.5, zorder=48) # adds coastline (since that isn't included)\n",
        "    ax.add_feature(cfeature.BORDERS, edgecolor='k', facecolor=None, linewidth=0.5, alpha=0.75, zorder=48) # adds country borders (since that isn't included)\n",
        "    cx.add_basemap(ax, source=cx.providers.CartoDB.DarkMatterOnlyLabels, zoom=zoom, zorder=49, attribution=False, crs=proj)\n",
        "    maplayertext1 = \"Map tiles: \" + str(cx.providers.CartoDB.DarkMatter.attribution)\n",
        "\n",
        "  elif map_theme == \"ESRI Light Grey\":\n",
        "    cx.add_basemap(ax, source='https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}', zoom=zoom, zorder=-1, crs=proj, attribution=False)\n",
        "    ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', edgecolor='#e4e4e4', facecolor='none', linewidth=0.5, zorder=48)) # adds state borders\n",
        "    ax.add_feature(cfeature.COASTLINE, edgecolor='#efefef', facecolor=None, linewidth=0.25, alpha=0.5, zorder=48) # adds coastline (since that isn't included)\n",
        "    ax.add_feature(cfeature.BORDERS, edgecolor='#efefef', facecolor=None, linewidth=0.5, alpha=0.75, zorder=48) # adds country borders (since that isn't included)\n",
        "    cx.add_basemap(ax, source='https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Reference/MapServer/tile/{z}/{y}/{x}', zoom=zoom, zorder=49, crs=proj, attribution=False)\n",
        "    maplayertext1 = \"Map tiles: Esri, HERE, Garmin, FAO, NOAA, USGS, © OpenStreetMap contributors, and the GIS User Community\"\n",
        "\n",
        "  elif map_theme == \"ESRI Dark Grey\":\n",
        "    cx.add_basemap(ax, source='https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Dark_Gray_Base/MapServer/tile/{z}/{y}/{x}', zorder=-1, crs=proj, attribution=False)\n",
        "    ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', edgecolor='#363638', facecolor='none', linewidth=0.5, zorder=48, alpha=0.5)) # adds state borders\n",
        "    ax.add_feature(cfeature.COASTLINE, edgecolor='#3f3f3f', facecolor=None, linewidth=0.25, alpha=0.5, zorder=48) # adds coastline (since that isn't included)\n",
        "    ax.add_feature(cfeature.BORDERS, edgecolor='#3f3f3f', facecolor=None, linewidth=0.5, alpha=0.75, zorder=48) # adds country borders (since that isn't included)\n",
        "    cx.add_basemap(ax, source='https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Dark_Gray_Reference/MapServer/tile/{z}/{y}/{x}', zoom=zoom, zorder=49, crs=proj, attribution=False)\n",
        "    maplayertext1 = \"Map tiles: Esri, HERE, Garmin, © OpenStreetMap contributors, and the GIS User Community\"\n",
        "\n",
        "  if cwa_opt:\n",
        "    try:\n",
        "      if os.path.exists(\"shp/w_22mr22.shp\"):\n",
        "        pass\n",
        "      else:\n",
        "        cwa_url = \"https://www.weather.gov/source/gis/Shapefiles/WSOM/w_22mr22.zip\"\n",
        "        cmd_mkshp = 'mkdir shp'\n",
        "        os.system(cmd_mkshp)\n",
        "        urlretrieve(cwa_url, \"shp/nws_cwa_outlines.zip\")\n",
        "        cmd_uz = 'unzip shp/nws_cwa_outlines.zip -d shp'\n",
        "        os.system(cmd_uz)\n",
        "      cwa_feature = ShapelyFeature(Reader(\"shp/w_22mr22.shp\").geometries(),ccrs.PlateCarree(), edgecolor='dimgrey', facecolor='none', linewidth=1.0, linestyle='dotted', zorder=3)\n",
        "      ax.add_feature(cwa_feature)\n",
        "    except:\n",
        "      print(\"    ❌ Aw shucks, no CWA boundaries for you. Sorry bout that.\")\n",
        "\n",
        "\n",
        "### Set up the plotting function (the thing that actually generates the finished figure)\n",
        "def drawmap(DATA,TITLESTRING,PROD,UNITS,LEVS):\n",
        "  global var\n",
        "  if not overview_opt:\n",
        "    print(\" > Finishing up map and adding legend\")\n",
        "  F = plt.gcf()  # Gets the current figure\n",
        "  ax = plt.gca()  # Gets the current axes\n",
        "  proxy = [mpatches.Patch(color = pc.get_facecolor()[0]) for pc in DATA.collections]\n",
        "  LLABS = []\n",
        "  for i in range(0, len(LEVS)-1):\n",
        "    label = str(LEVS[i])+\"-\"+str(LEVS[i+1])+UNITS\n",
        "    LLABS.append(label)\n",
        "\n",
        "  if ob_opt:\n",
        "    if nbm_var == \"PMxMnT\":\n",
        "      if plot_ob_misses:\n",
        "        LLABS.append(icon_dict[element][\"Ob_miss\"][\"label\"])\n",
        "        proxy.append(plt.scatter([], [], marker=\"o\", color=\"#000000\", edgecolor=icon_dict[element][\"Ob_miss\"][\"color\"], linewidth=0.5, s=20, transform=proj))\n",
        "      LLABS.append(icon_dict[element][\"Ob_hit\"][\"label\"])\n",
        "      proxy.append(plt.scatter([], [], marker=\"o\", color=icon_dict[element][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.5, s=20, transform=proj))\n",
        "    elif nbm_var == \"PQPF\":\n",
        "      LLABS.append(icon_dict[ptype][\"Ob\"][\"label\"])\n",
        "      proxy.append(plt.scatter([], [], marker=get_marker(icon_dict[ptype][\"Ob\"][\"icon\"]), color=icon_dict[ptype][\"Ob\"][\"color\"], edgecolor='w', linewidth=0.1, s=35,transform=proj))\n",
        "    elif nbm_var == \"PMxWind\":\n",
        "      if plot_ob_misses:\n",
        "        LLABS.append(icon_dict[nbm_var][\"Ob_miss\"][\"label\"])\n",
        "        proxy.append(plt.scatter([], [], marker=\"o\", color=\"#000000\", edgecolor=icon_dict[nbm_var][\"Ob_miss\"][\"color\"], linewidth=0.5, s=20, transform=proj))\n",
        "      LLABS.append(icon_dict[nbm_var][\"Ob_hit\"][\"label\"])\n",
        "      proxy.append(plt.scatter([], [], marker=\"o\", color=icon_dict[nbm_var][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.5, s=20, transform=proj))\n",
        "\n",
        "  if lsr_opt:\n",
        "    if nbm_var == \"PQPF\":\n",
        "      if ptype == \"qpf\":\n",
        "        LLABS.append(icon_dict[ptype][\"LSR\"][\"heavyrain\"][\"label\"])\n",
        "        proxy.append(plt.scatter([], [], marker=get_marker(icon_dict[ptype][\"LSR\"][\"heavyrain\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"heavyrain\"][\"color\"], edgecolor='w', linewidth=0.1, s=35,transform=proj))\n",
        "        LLABS.append(icon_dict[ptype][\"LSR\"][\"flooding\"][\"label\"])\n",
        "        proxy.append(plt.scatter([], [], marker=get_marker(icon_dict[ptype][\"LSR\"][\"flooding\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"flooding\"][\"color\"], edgecolor='w', linewidth=0.1, s=35,transform=proj))\n",
        "      else:\n",
        "        LLABS.append(icon_dict[ptype][\"LSR\"][\"label\"])\n",
        "        proxy.append(plt.scatter([], [], marker=get_marker(icon_dict[ptype][\"LSR\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"color\"], edgecolor='w', linewidth=0.1, s=35, transform=proj))\n",
        "    elif nbm_var == \"PMxWind\":\n",
        "      LLABS.append(icon_dict[nbm_var][\"LSR\"][\"label\"])\n",
        "      proxy.append(plt.scatter([], [], marker=get_marker(icon_dict[nbm_var][\"LSR\"][\"icon\"]), color=icon_dict[nbm_var][\"LSR\"][\"color\"], edgecolor='w', linewidth=0.1, s=35, transform=proj))\n",
        "\n",
        "  proxy = proxy[::-1]\n",
        "  LLABS = LLABS[::-1]\n",
        "\n",
        "  l = plt.legend(handles=proxy, labels=LLABS, fontsize='5',fancybox=True, loc=LLOC)\n",
        "  for text in l.get_texts():\n",
        "    text.set_color(font_color)\n",
        "  lframe = l.get_frame()\n",
        "  lframe.set_color(face_color)\n",
        "  l.set_zorder(100)\n",
        "\n",
        "  #artists, labels = DATA.legend_elements()\n",
        "\n",
        "  nbm_init_title = nbm_init.strftime('%HZ %d-%b-%Y')\n",
        "  nbm_text = r\"$\\bf{\"+'NBM'+\"}$\" +' (National Blend of Models) · ' +r\"$\\bf{\"+\"Init\"+\"}$\" +\" \"+ nbm_init_title\n",
        "  if nbm_var == \"PMxMnT\" or nbm_var == \"PMxWind\":\n",
        "    #nbm_text = f'NBM · National Blend of Models · Init {nbm_init_title}'\n",
        "    valid_date_title = valid_date_start.strftime('%A\\n %b %d %Y')\n",
        "  else:\n",
        "    valid_date_title = valid_date.strftime('%HZ %a\\n %b %d %Y')\n",
        "\n",
        "  plt.text(0.0,1.03, TITLESTRING, transform=ax.transAxes, ha='left', va='bottom', weight='bold', fontsize=9, color=font_color)\n",
        "  plt.text(0.0,1.03, nbm_text, transform=ax.transAxes, ha='left', va='top', fontsize=6, color='grey')\n",
        "  plt.text(1.0,1.03,valid_date_title, transform=ax.transAxes, ha='right', va='center', weight='bold', fontsize=7, color=font_color)\n",
        "\n",
        "  if maplayertext1 != \"\":\n",
        "    if maplayertext2 !=\"\":\n",
        "        plt.text(0.5, -0.009, '%s, %s // Made with Wx4Colab' % (maplayertext1,maplayertext2),\n",
        "                transform = ax.transAxes, ha='center', va ='top',fontsize=4,color=font_color, style='italic', zorder=99,\n",
        "                bbox=dict(facecolor=face_color,edgecolor='none', pad=1.8, alpha=0.65))\n",
        "    else:\n",
        "      plt.text(0.5, -0.009, '%s // Made with Wx4Colab' % (maplayertext1),\n",
        "                transform = ax.transAxes, ha='center', va ='top',fontsize=4,color=font_color, style='italic', zorder=99,\n",
        "                bbox=dict(facecolor=face_color,edgecolor='none', pad=1.8, alpha=0.65))\n",
        "  if var == \"percentile\":\n",
        "    pthresh = threshold_percentile\n",
        "  else:\n",
        "    if nbm_var == \"PQPF\":\n",
        "      if ptype==\"qpf\":\n",
        "        pthresh = threshold_qpf\n",
        "      elif ptype==\"snow\":\n",
        "        pthresh = threshold_snow\n",
        "      elif ptype==\"ice\":\n",
        "        pthresh = threshold_ice\n",
        "    elif nbm_var == \"PMxWind\":\n",
        "      pthresh = wind_threshold\n",
        "    elif nbm_var == \"PMxMnT\":\n",
        "      pthresh = threshold_temp\n",
        "\n",
        "  if nbm_var == \"PMxMnT\":\n",
        "    desc = temp\n",
        "    var = \"prob\"\n",
        "    flengthtitle = \"18\"\n",
        "  elif nbm_var == \"PQPF\":\n",
        "    desc = ptype\n",
        "    flengthtitle = str(forecast_length)\n",
        "  elif nbm_var == \"PMxWind\":\n",
        "    desc = f\"max{wind_type}\"\n",
        "    flengthtitle = \"24\"\n",
        "  if dom==\"Custom\" and custom_name:\n",
        "    custom_name_sanitized = custom_name.replace(\" \",\"\")\n",
        "    file_id = '%s_%s%s%s_%s_%sT%s' % (custom_name_sanitized, var, desc, flengthtitle, pthresh, nbm_init.strftime('%Y%m%d%H'), valid_date.strftime('%Y%m%d%H'))\n",
        "  else:\n",
        "    file_id = '%s_%s%s%s_%s_%sT%s' % (dom, var, desc, flengthtitle, pthresh, nbm_init.strftime('%Y%m%d%H'), valid_date.strftime('%Y%m%d%H'))\n",
        "\n",
        "  filename = '%s.png' % (file_id)\n",
        "  print(f' > Saving plot as {filename}')\n",
        "  plt.savefig(filename,bbox_inches='tight', facecolor=face_color) # Saves the figure with small margins\n",
        "\n",
        "  plt.show()\n",
        "\n",
        "\n",
        "\n",
        "# These set up how the plot looks like based on the variable you have chosen to plot.\n",
        "##### Snow Plots #######################################################################################################\n",
        "def plot_snow_prob():\n",
        "  print(\"Making snow probability plot\")\n",
        "  blankmap()\n",
        "  probsnow = nbm\n",
        "  units = r\"%\"\n",
        "  BuPu_Stacked =[\"#E1E9EF\", \"#C3D6E9\", \"#A4C5E3\", \"#86B4DD\", \"#66A1D7\",\n",
        "                  \"#CA56A6\", \"#B54BA3\", \"#9E42A0\", \"#88389E\", \"#6C2B9B\"]\n",
        "  BuPu_Converging = [\"#CBDDF0BF\",\"#8DC1DDBF\",\"#509BCCBF\",\"#2273B6BF\",\"#074991BF\",\n",
        "                      \"#3F007DBF\",\"#62419CBF\",\"#7E79B9BF\",\"#A29ECBBF\",\"#C9CAE2BF\"]\n",
        "  BuPu_Continuous = [\"#d4ecf9BF\",\"#add1ebBF\",\"#85b8ddBF\",\"#679ccdBF\",\"#4a83bfBF\",\n",
        "                \"#376eb1BF\",\"#3d61a8BF\",\"#42519dBF\",\"#3d3993BF\",\"#321486BF\"]\n",
        "  GyBuPk = [\"#5e5e5e\", \"#5a6375\", \"#56678c\", \"#536b9f\", \"#506fb5\",\n",
        "            \"#5b75c1\", \"#757cc7\", \"#9586d0\", \"#b48fd7\", \"#d99ae0\"]\n",
        "  if \"Dark\" in map_theme:\n",
        "    cmap = GyBuPk\n",
        "  else:\n",
        "    cmap = BuPu_Stacked\n",
        "  CLEVS = [5,10,20,30,40,50,60,70,80,90,100]\n",
        "  print(\" > Contouring NBM\")\n",
        "  probsnow_plot = plt.contourf(nbmlons,nbmlats,probsnow,CLEVS,colors=cmap,extend='neither', alpha=0.75, transform=latloncrs, antialiased = True, transform_first=True)\n",
        "  #probsnow_plot = plt.contourf(x,y,probsnow,CLEVS,cmap=plt.get_cmap('BuPu'),extend='neither',alpha=0.75, transform=datacrs, antialiased = True)\n",
        "  #uncomment the following lines if you want contours and labels\n",
        "  #contours = plt.contour(x,y,probsnow,CLEVS,colors='silver',linewidths=1., alpha=0.35, transform=latloncrs)\n",
        "  #clabs = plt.clabel(contours, inline=1, fontsize=8, fmt='%1.0f', use_clabeltext=True, colors=font_color)\n",
        "  #plt.setp(clabs, fontweight='bold', path_effects=[PathEffects.withStroke(linewidth=1.4,foreground=face_color)])\n",
        "\n",
        "  if lsr_opt:\n",
        "    print(\" > Getting LSRs\")\n",
        "    #lsr_lats, lsr_lons, lsr_mag = [], [], []\n",
        "    lsr_start = valid_date - timedelta (hours=int(forecast_length))\n",
        "    lsr_start = lsr_start.strftime('%Y-%m-%dT%H:%M')\n",
        "    lsr_end = valid_date.strftime('%Y-%m-%dT%H:%M')\n",
        "    lsr_url = \"https://mesonet.agron.iastate.edu/cgi-bin/request/gis/lsr.py?wfo[]=ALL&sts=\"+lsr_start+\"Z&ets=\"+lsr_end+\"Z&fmt=csv\"\n",
        "    csv_name = \"lsr_\"+lsr_start+\"_\"+lsr_end+\".csv\"\n",
        "\n",
        "    try:\n",
        "      file_exists = os.path.exists(csv_name)\n",
        "      if file_exists:\n",
        "        pass\n",
        "      else:\n",
        "        urlretrieve(lsr_url, csv_name)\n",
        "\n",
        "      LSRs = pd.read_csv(csv_name, usecols= ['LAT','LON','MAG','TYPECODE','TYPETEXT'], on_bad_lines='skip')\n",
        "      LSRs['MAG'] = LSRs.MAG.replace('None','0').astype(float)\n",
        "      LSRs['LAT'] = LSRs.LAT.astype(float)\n",
        "      LSRs['LON'] = LSRs.LON.astype(float)\n",
        "      #LSRs[(LSRs.TYPECODE=='5')&(LSRs.MAG>=0.1)].plot.scatter(x='LON', y='LAT', c='MAG', vmin=0, vmax=1, ec='red', s=100, transform=proj)\n",
        "      LSRs_filtered = LSRs[(LSRs.TYPECODE == 'S') & (LSRs.MAG >= float(threshold_snow))]\n",
        "      lsr_lats= LSRs_filtered['LAT'].to_numpy()\n",
        "      lsr_lons= LSRs_filtered['LON'].to_numpy()\n",
        "      lsr_mag= LSRs_filtered['MAG'].to_numpy()\n",
        "      plt.scatter(lsr_lons, lsr_lats, marker=get_marker(icon_dict[ptype][\"LSR\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"color\"], edgecolor='w', linewidth=0.1, s=15,transform=latloncrs, zorder=3)\n",
        "      print(\"   ✅ ...and plotted!\")\n",
        "\n",
        "    except:\n",
        "      print(\"   ❌ No LSRs for you\")\n",
        "\n",
        "  title = 'Chance of More Than '+str(threshold_snow)+'\\\" of Snow in '+forecast_length+'hrs'\n",
        "  prodid = 'probsnow'+forecast_length+'_'+str(threshold_snow)\n",
        "  drawmap(probsnow_plot, title, prodid, units, CLEVS)\n",
        "\n",
        "def plot_snow_percentile():\n",
        "  print(\"Making snow percentile plot\")\n",
        "  blankmap()\n",
        "  snowp = nbm\n",
        "  SNOW_COLS=[\n",
        "      \"#bdd7e7\", # 0.1\" - 1\"\n",
        "      \"#6bafd6\", # 1\" - 2\"\n",
        "      \"#3182bd\", # 2\" - 3\"\n",
        "      \"#08529c\", # 3\" - 4\"\n",
        "      \"#082694\", # 4\" - 6\"\n",
        "      \"#ffff96\", # 6\" - 8\"\n",
        "      \"#ffc400\", # 8\" - 12\"\n",
        "      \"#ff8800\", # 12\" - 18\"\n",
        "      \"#db1200\", # 18\" - 24\"\n",
        "      \"#9e0000\", # 24\" - 30\"\n",
        "      \"#690000\", # 30\" - 36\"\n",
        "      \"#360000\", # 36\" - 48\"\n",
        "  ]\n",
        "  #        \"#ccccff\", # 48\" - 60\"\n",
        "  #        \"#9f8cd8\", # 60\" - 72\"\n",
        "  #        \"#7c52a5\", # 72\" - 96\"\n",
        "  #        \"#551c72\", # 96\" - 120\"\n",
        "  #        \"#2e003d\", # 120+\"\n",
        "  #    ]\n",
        "  SNOW_LEVS = [0.1,1,2,3,4,6,8,12,18,24,30,36,48]\n",
        "  #SNOW_LEVS = [0.1,1,2,3,4,6,8,12,18,24,30,36,48,60,72,96,120]\n",
        "  print(\" > Contouring NBM\")\n",
        "  SNOW = plt.contourf(nbmlons,nbmlats,snowp,SNOW_LEVS,colors=SNOW_COLS,extend='neither', alpha=0.75, transform=latloncrs, antialiased = True, transform_first=True)\n",
        "  units = \" in\"\n",
        "  title = str(threshold_percentile)+perc_suffix_dict[threshold_percentile]+' Percentile '+forecast_length+'hr Snow'\n",
        "  prodid = 'snowp'+forecast_length+'_'+str(threshold_snow)\n",
        "\n",
        "  drawmap(SNOW, title, prodid, units, SNOW_LEVS)\n",
        "\n",
        "\n",
        "\n",
        "###### PQPF plots ######################################################################################################\n",
        "\n",
        "def plot_qpf_prob():\n",
        "  print(\"Making QPF probability plot\")\n",
        "  blankmap()\n",
        "  probqpf = nbm\n",
        "  units = r\"%\"\n",
        "  GuBu_stacked = [\"#F5F7D4\",\"#DCE9A7\",\"#B9DB88\",\"#92C86A\",\"#6DB54C\",\n",
        "                \"#44A236\",\"#269130\",\"#1C753A\",\"#1C5E44\",\"#254944\"]\n",
        "  GyGrYGn = [\"#5e5e5e\", \"#5b6853\", \"#586f4b\", \"#567642\", \"#557e39\",\n",
        "            \"#5e8a2f\", \"#7a9e23\", \"#98af3c\", \"#b2c440\", \"#ccd945\"]\n",
        "  if \"Dark\" in map_theme:\n",
        "    cmap = GyGrYGn\n",
        "  else:\n",
        "    cmap = GuBu_stacked\n",
        "  CLEVS = [5,10,20,30,40,50,60,70,80,90,100]\n",
        "  print(\" > Contouring NBM\")\n",
        "  probqpf_plot = plt.contourf(nbmlons,nbmlats,probqpf,CLEVS,colors=cmap,extend='neither', alpha=0.75, transform=latloncrs, antialiased = True, transform_first=True)\n",
        "  #uncomment the following lines if you want contours and labels\n",
        "  #contours = plt.contour(x,y,probqpf,CLEVS,colors='silver',linewidths=1., alpha=0.35, transform=latloncrs)\n",
        "  #clabs = plt.clabel(contours, inline=1, fontsize=8, fmt='%1.0f', use_clabeltext=True, colors=font_color)\n",
        "  #plt.setp(clabs, fontweight='bold', path_effects=[PathEffects.withStroke(linewidth=1.4,foreground=face_color)])\n",
        "\n",
        "  if ob_opt:\n",
        "    print(\" > Getting obs\")\n",
        "    bbox_str = str(west)+\",\"+str(south)+\",\"+str(east)+\",\"+str(north)\n",
        "    obs_start = valid_date - timedelta (hours=int(forecast_length))\n",
        "    obs_start = obs_start.strftime('%Y%m%d%H%M')\n",
        "    obs_end = valid_date.strftime('%Y%m%d%H%M')\n",
        "    obs_url = \"https://api.synopticdata.com/v2/stations/precipitation?&token=\"+synoptic_token+\"&start=\"+obs_start+\"&end=\"+\\\n",
        "              obs_end+\"&pmode=totals&obtimezone=utc&units=precip|in&\"+network_string+\"&bbox=\"+bbox_str+\"&fields=latitude,longitude\"\n",
        "    json_name = \"obs_\"+dom+ \"_\"+obs_start+\"_\"+obs_end+\"_\"+str(south).replace('.','')+\".json\"\n",
        "    obs_lats = []\n",
        "    obs_lons = []\n",
        "    obs_stid = []\n",
        "    obs_mag = []\n",
        "    try:\n",
        "      json_exists = os.path.exists(json_name)\n",
        "      if json_exists:\n",
        "        pass\n",
        "      else:\n",
        "        urlretrieve(obs_url, json_name)\n",
        "      with open(json_name) as json_file:\n",
        "          obs_json = json.load(json_file)\n",
        "          for stn in obs_json[\"STATION\"]:\n",
        "            olat = stn[\"LATITUDE\"]\n",
        "            olon = stn[\"LONGITUDE\"]\n",
        "            #ptotal = None\n",
        "            if stn[\"STID\"] is None:\n",
        "              stid = \"N0N3\"\n",
        "            else:\n",
        "              stid = stn[\"STID\"]\n",
        "\n",
        "            if \"precipitation\" in stn[\"OBSERVATIONS\"]:\n",
        "\n",
        "              if \"total\" in stn[\"OBSERVATIONS\"][\"precipitation\"][0]:\n",
        "                ptotal = stn[\"OBSERVATIONS\"][\"precipitation\"][0][\"total\"]\n",
        "\n",
        "                if ptotal >= (float(threshold_qpf)):\n",
        "                  obs_lats.append(float(olat))\n",
        "                  obs_lons.append(float(olon))\n",
        "                  obs_stid.append(stid)\n",
        "                  obs_mag.append(ptotal)\n",
        "\n",
        "                #obs_lats.append(float(olat))\n",
        "                #obs_lons.append(float(olon))\n",
        "                #obs_stid.append(stid)\n",
        "                #obs_mag.append(ptotal)\n",
        "\n",
        "      #obs_df = pd.DataFrame()\n",
        "      #obs_df[\"lat\"] = obs_lats\n",
        "      #obs_df[\"lon\"] = obs_lons\n",
        "      #obs_df[\"stid\"] = obs_stid\n",
        "      #obs_df[\"ptotal\"] = obs_mag\n",
        "      #csv_filename = f\"obs_precip_{obs_end}.csv\"\n",
        "      #obs_df.to_csv(csv_filename)\n",
        "      plt.scatter(obs_lons, obs_lats, marker=get_marker(icon_dict[ptype][\"Ob\"][\"icon\"]), color=icon_dict[ptype][\"Ob\"][\"color\"], edgecolor='w', linewidth=0.1, s=15,transform=latloncrs, zorder=3)\n",
        "      print(\"   ✅ ...and plotted!\")\n",
        "    except:\n",
        "      (\"   ❌ No obs for you\")\n",
        "\n",
        "  if lsr_opt:\n",
        "    print(\" > Getting LSRs\")\n",
        "    #lsr_lats, lsr_lons, lsr_mag = [], [], []\n",
        "    lsr_start = valid_date - timedelta (hours=int(forecast_length))\n",
        "    lsr_start = lsr_start.strftime('%Y-%m-%dT%H:%M')\n",
        "    lsr_end = valid_date.strftime('%Y-%m-%dT%H:%M')\n",
        "    lsr_url = \"https://mesonet.agron.iastate.edu/cgi-bin/request/gis/lsr.py?wfo[]=ALL&sts=\"+lsr_start+\"Z&ets=\"+lsr_end+\"Z&fmt=csv\"\n",
        "    csv_name = \"lsr_\"+lsr_start+\"_\"+lsr_end+\".csv\"\n",
        "\n",
        "    try:\n",
        "      file_exists = os.path.exists(csv_name)\n",
        "      if file_exists:\n",
        "        pass\n",
        "      else:\n",
        "        urlretrieve(lsr_url, csv_name)\n",
        "      LSRs = pd.read_csv(csv_name, usecols= ['LAT','LON','MAG','TYPECODE','TYPETEXT'], on_bad_lines='skip')\n",
        "      LSRs['MAG'] = LSRs.MAG.replace('None','0').astype(float)\n",
        "      LSRs['LAT'] = LSRs.LAT.astype(float)\n",
        "      LSRs['LON'] = LSRs.LON.astype(float)\n",
        "      LSRs_HeavyRain = LSRs[(LSRs.TYPECODE == 'R') & (LSRs.MAG >= float(threshold_qpf))]\n",
        "      heavyrain_lats= LSRs_HeavyRain['LAT'].to_numpy()\n",
        "      heavyrain_lons= LSRs_HeavyRain['LON'].to_numpy()\n",
        "      heavyrain_mag= LSRs_HeavyRain['MAG'].to_numpy()\n",
        "      plt.scatter(heavyrain_lons, heavyrain_lats, marker=get_marker(icon_dict[ptype][\"LSR\"][\"heavyrain\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"heavyrain\"][\"color\"], edgecolor=None, linewidth=0.0, s=35,transform=latloncrs, zorder=4)\n",
        "\n",
        "      LSRs_flooding = LSRs[(LSRs.TYPECODE == 'E') | (LSRs.TYPECODE == 'F')]\n",
        "      flood_lats= LSRs_flooding['LAT'].to_numpy()\n",
        "      flood_lons= LSRs_flooding['LON'].to_numpy()\n",
        "      plt.scatter(flood_lons, flood_lats, marker=get_marker(icon_dict[ptype][\"LSR\"][\"flooding\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"flooding\"][\"color\"], edgecolor=None, linewidth=0.0, s=50,transform=latloncrs, zorder=5)\n",
        "      print(\"   ✅ ...and plotted!\")\n",
        "\n",
        "\n",
        "    except:\n",
        "      (\"   ❌ No LSRs for you\")\n",
        "\n",
        "\n",
        "  title = 'Chance of More Than '+str(threshold_qpf)+'\\\" of QPF in '+forecast_length+'hrs'\n",
        "  prodid = 'probqpf'+forecast_length+'_'+str(threshold_qpf)\n",
        "\n",
        "  drawmap(probqpf_plot, title, prodid, units, CLEVS)\n",
        "\n",
        "def plot_qpf_percentile():\n",
        "  print(\"Making QPF percentile plot\")\n",
        "  blankmap()\n",
        "  qpf = nbm\n",
        "  QPF_COLS = [\"#c7e9c0\",\"#a1d99b\", \"#74c476\", \"#31a353\", \"#006d2c\", \"#fff98a\", \"#ffcc4f\", \"#fe8d3c\", \\\n",
        "    \"#fc4e2a\", \"#d61a1d\", \"#ad0025\", \"#700025\", \"#3b0030\", \"#4d0073\", \"#ffdbff\"]\n",
        "  QPF_LEVS = [0.01,0.10,0.25,0.50,0.75,1.0,1.5,2.0,3.0,4.0,6.0,8.0,10.0,15.0,20.0,30.0]\n",
        "  print(\" > Contouring NBM\")\n",
        "  PQPF = plt.contourf(nbmlons,nbmlats,qpf,QPF_LEVS,colors=QPF_COLS,extend='neither', alpha=0.75, transform=latloncrs,antialiased = True, transform_first=True)\n",
        "  units = \" in\"\n",
        "  title = str(threshold_percentile)+perc_suffix_dict[threshold_percentile]+' Percentile '+forecast_length+'hr QPF'\n",
        "  prodid = 'pqpf'+forecast_length+'_'+str(threshold_percentile)\n",
        "\n",
        "  drawmap(PQPF, title, prodid, units, QPF_LEVS)\n",
        "\n",
        "\n",
        "def plot_qpf_percentile_overview():\n",
        "  with io.capture_output() as captured:\n",
        "    blankmap()\n",
        "  qpf = nbm\n",
        "  QPF_COLS = [\"#c7e9c0\",\"#a1d99b\", \"#74c476\", \"#31a353\", \"#006d2c\", \"#fff98a\", \"#ffcc4f\", \"#fe8d3c\", \\\n",
        "    \"#fc4e2a\", \"#d61a1d\", \"#ad0025\", \"#700025\", \"#3b0030\", \"#4d0073\", \"#ffdbff\"]\n",
        "  QPF_LEVS = [0.01,0.10,0.25,0.50,0.75,1.0,1.5,2.0,3.0,4.0,6.0,8.0,10.0,15.0,20.0,30.0]\n",
        "  PQPF = plt.contourf(nbmlons,nbmlats,qpf,QPF_LEVS,colors=QPF_COLS,extend='neither', alpha=0.75, transform=latloncrs,antialiased = True, transform_first=True)\n",
        "  units = \" in\"\n",
        "  title = str(threshold_percentile)+perc_suffix_dict[threshold_percentile]+' Percentile '+forecast_length+'hr QPF'\n",
        "  prodid = 'pqpf'+forecast_length+'_'+str(threshold_percentile)\n",
        "\n",
        "  drawmap(PQPF, title, prodid, units, QPF_LEVS)\n",
        "\n",
        "\n",
        "\n",
        "### PICE plots #########################################################################################################\n",
        "\n",
        "def plot_ice_prob():\n",
        "  print(\"Making ice probability plot\")\n",
        "  blankmap()\n",
        "  probice= nbm\n",
        "  units = r\"%\"\n",
        "  #prob_cols = [\"#E4F0F9\",\"#B7D4EA\",\"#69ADD5\",\"#2E7EBC\",\"#074991\",\n",
        "  #             \"#F994B2\",\"#F15A9E\",\"#CD238E\",\"#99017B\",\"#56006D\"]\n",
        "  GyMgYl = [\"#5e5e5e\",\"#745c77\",\"#855b8a\",\"#975b9f\",\"#a95ab3\",\n",
        "            \"#b56cb4\",\"#bf89ab\",\"#cbaaa2\",\"#d7c99b\",\"#e4f287\"]\n",
        "  YlPkMg = [\"#e4f287\",\"#ded967\",\"#d8c09a\",\"#d2aa9e\",\"#cd95a2\",\n",
        "            \"#c17f9f\",\"#b06b97\",\"#9c558e\",\"#894085\",\"#74287c\"]\n",
        "  if \"Dark\" in map_theme:\n",
        "    cmap = GyMgYl\n",
        "  else:\n",
        "    cmap = YlPkMg\n",
        "  CLEVS = [5,10,20,30,40,50,60,70,80,90,100]\n",
        "  print(\" > Contouring NBM\")\n",
        "  probice_plot = plt.contourf(nbmlons,nbmlats,probice,CLEVS,colors=cmap, extend='neither', alpha=0.75, transform=latloncrs, antialiased = True, transform_first=True)\n",
        "  #uncomment the following lines if you want contours and labels\n",
        "  #contours = plt.contour(x,y,probice,CLEVS,colors='silver',linewidths=1., alpha=0.35, transform=latloncrs)\n",
        "  #clabs = plt.clabel(contours, inline=1, fontsize=8, fmt='%1.0f', use_clabeltext=True, colors=font_color)\n",
        "  #plt.setp(clabs, fontweight='bold', path_effects=[PathEffects.withStroke(linewidth=1.4,foreground=face_color)])\n",
        "\n",
        "  if lsr_opt:\n",
        "    print(\" > Getting LSRs\")\n",
        "    lsr_start = valid_date - timedelta (hours=int(forecast_length))\n",
        "    lsr_start = lsr_start.strftime('%Y-%m-%dT%H:%M')\n",
        "    lsr_end = valid_date.strftime('%Y-%m-%dT%H:%M')\n",
        "    lsr_url = \"https://mesonet.agron.iastate.edu/cgi-bin/request/gis/lsr.py?wfo[]=ALL&sts=\"+lsr_start+\"Z&ets=\"+lsr_end+\"Z&fmt=csv\"\n",
        "    csv_name = \"lsr_\"+lsr_start+\"_\"+lsr_end+\".csv\"\n",
        "\n",
        "    try:\n",
        "      file_exists = os.path.exists(csv_name)\n",
        "      if file_exists:\n",
        "        pass\n",
        "      else:\n",
        "        urlretrieve(lsr_url, csv_name)\n",
        "\n",
        "      LSRs = pd.read_csv(csv_name, usecols= ['LAT','LON','MAG','TYPECODE','TYPETEXT'], on_bad_lines='skip')\n",
        "      LSRs['MAG'] = LSRs.MAG.replace('None','0').astype(float)\n",
        "      LSRs['LAT'] = LSRs.LAT.astype(float)\n",
        "      LSRs['LON'] = LSRs.LON.astype(float)\n",
        "      LSRs_filtered = LSRs[(LSRs.TYPECODE == '5') & (LSRs.MAG >= float(threshold_ice))]\n",
        "      lsr_lats= LSRs_filtered['LAT'].to_numpy()\n",
        "      lsr_lons= LSRs_filtered['LON'].to_numpy()\n",
        "      lsr_mag= LSRs_filtered['MAG'].to_numpy()\n",
        "      plt.scatter(lsr_lons, lsr_lats, marker=get_marker(icon_dict[ptype][\"LSR\"][\"icon\"]), color=icon_dict[ptype][\"LSR\"][\"color\"], edgecolor=None, linewidth=0.0, s=15,transform=latloncrs, zorder=3)\n",
        "      print(\"   ✅ ...and plotted!\")\n",
        "    except:\n",
        "      (\"   ❌ No LSRs for you\")\n",
        "\n",
        "\n",
        "\n",
        "  title = 'Chance of More Than '+str(threshold_ice)+' Inches of Ice in '+forecast_length+'hrs'\n",
        "  prodid = 'probice'+forecast_length+'_'+str(threshold_ice)\n",
        "\n",
        "  drawmap(probice_plot, title, prodid, units, CLEVS)\n",
        "\n",
        "def plot_ice_percentile():\n",
        "  print(\"Making ice percentile plot\")\n",
        "  blankmap()\n",
        "  ice = nbm\n",
        "  #QPF_COLS = [\"#c7e9c0\",\"#a1d99b\", \"#74c476\", \"#31a353\", \"#006d2c\", \"#fff98a\", \"#ffcc4f\", \"#fe8d3c\", \\\n",
        "  #  \"#fc4e2a\", \"#d61a1d\", \"#ad0025\", \"#700025\", \"#3b0030\", \"#4d0073\", \"#ffdbff\"]\n",
        "  ICE_LEVS = [0.01,0.10,0.22,0.30,0.40,0.50,0.75,1.0,1.25,1.5,2.0]\n",
        "  print(\" > Contouring NBM\")\n",
        "  PICE = plt.contourf(nbmlons,nbmlats,ice,ICE_LEVS,cmap=plt.get_cmap('spring_r'),extend='neither', alpha=0.75, transform=latloncrs,antialiased = True, transform_first=True)\n",
        "  #uncomment the following lines if you want contours and labels\n",
        "  #contours = plt.contour(x,y,ice,ICE_LEVS,colors='silver',linewidths=1., alpha=0.35, transform=latloncrs)\n",
        "  #clabs = plt.clabel(contours, inline=1, fontsize=8, fmt='%1.0f', use_clabeltext=True, colors=font_color)\n",
        "  #plt.setp(clabs, fontweight='bold', path_effects=[PathEffects.withStroke(linewidth=1.4,foreground=face_color)])\n",
        "  units = \" in\"\n",
        "  title = str(threshold_percentile)+perc_suffix_dict[threshold_percentile]+' Percentile '+forecast_length+'hr Ice'\n",
        "  prodid = 'pice'+forecast_length+'_'+str(threshold_ice)\n",
        "\n",
        "  drawmap(PICE, title, prodid, units, ICE_LEVS)\n",
        "\n",
        "\n",
        "### MaxT Plots ###\n",
        "def plot_maxt_prob():\n",
        "    print(\"Making prob MaxT plot...\")\n",
        "    blankmap()\n",
        "    print(\" > Contouring NBM\")\n",
        "    #probmaxt = nbm\n",
        "    #probmaxt[probmaxt < 5.] = -999\n",
        "    units = r\"%\"\n",
        "    CLEVS = [5,10,20,30,40,50,60,70,80,90,100]\n",
        "    GyBuWh = [\"#4D4D4D\", \"#4D566C\", \"#4E5E84\", \"#4E669A\", \"#4E6DB4\",\n",
        "            \"#5D7CC4\", \"#7D99D0\", \"#9BB2DB\", \"#BCCCE7\", \"#E0EAF5\"]\n",
        "    GyRdYl = [\"#4D4D4D\", \"#724747\", \"#914242\", \"#ad3e3e\", \"#c93939\",\n",
        "            \"#c94b46\", \"#e26e5f\", \"#ec9a7f\", \"#f6c59f\", \"#ffeebd\"]\n",
        "    BuPu_Continuous = [\"#d4ecf9BF\",\"#add1ebBF\",\"#85b8ddBF\",\"#679ccdBF\",\"#4a83bfBF\",\n",
        "                \"#376eb1BF\",\"#3d61a8BF\",\"#42519dBF\",\"#3d3993BF\",\"#321486BF\"]\n",
        "    if \"Cold\" in element:\n",
        "      if \"Dark\" in map_theme:\n",
        "        cmap = GyBuWh\n",
        "      else:\n",
        "        cmap_sel = plt.get_cmap('BuPu')\n",
        "        cmap=[cmap_sel(0.0), cmap_sel(0.1), cmap_sel(0.2), cmap_sel(0.3), cmap_sel(0.4), cmap_sel(0.5), cmap_sel(0.6), cmap_sel(0.7), cmap_sel(0.8), cmap_sel(0.9), cmap_sel(1.0)]\n",
        "    elif \"Hot\" in element:\n",
        "      if \"Dark\" in map_theme:\n",
        "        cmap = GyRdYl\n",
        "      else:\n",
        "        cmap_sel = plt.get_cmap('YlOrRd')\n",
        "        cmap=[cmap_sel(0.0), cmap_sel(0.1), cmap_sel(0.2), cmap_sel(0.3), cmap_sel(0.4), cmap_sel(0.5), cmap_sel(0.6), cmap_sel(0.7), cmap_sel(0.8), cmap_sel(0.9), cmap_sel(1.0)]\n",
        "    #plt.contourf(x, y, nbm, CLEVS, cmap=plt.get_cmap(colormap_temp), extend='neither', transform=datacrs, antialiased=True)\n",
        "    #nbmnull = nbm - nbm\n",
        "    probmaxt_plot = plt.contourf(nbmlons, nbmlats, nbm, CLEVS, colors=cmap, extend='neither',\n",
        "                                 alpha=0.75, transform=latloncrs, transform_first=True, antialiased=True)\n",
        "    #uncomment the following lines if you want contours and labels\n",
        "    #contours = plt.contour(x,y,nbm,CLEVS,colors='silver',linewidths=1., alpha=0.35, transform=latloncrs)\n",
        "    #clabs = plt.clabel(contours, inline=1, fontsize=8, fmt='%1.0f', use_clabeltext=True, colors=font_color)\n",
        "    #plt.setp(clabs, fontweight='bold', path_effects=[PathEffects.withStroke(linewidth=1.4,foreground=face_color)])\n",
        "\n",
        "\n",
        "    if ob_opt:\n",
        "      print(\" > Now looking for obs\")\n",
        "      bbox_str = str(west)+\",\"+str(south)+\",\"+str(east)+\",\"+str(north)\n",
        "      obs_url = \"https://api.synopticlabs.org/v2/stations/statistics?token=\"+synoptic_token+\"&vars=air_temp&start=\"+ \\\n",
        "                valid_date_start.strftime('%Y%m%d')+obs_start_hour+\"&end=\"+valid_date_end.strftime('%Y%m%d')+ \\\n",
        "                obs_end_hour+\"&units=temp%7Cf&within=1440&type=\"+ob_stat+\"&status=active\"+network_string+\"&bbox=\"+bbox_str+\"&fields=latitude,longitude\"\n",
        "      json_name = \"obs/Obs_MaxT_\"+valid_date_start.strftime('%Y%m%d')+obs_start_hour+\"_\"+valid_date_end.strftime('%Y%m%d')+obs_end_hour+\"_\"+dom+\".json\"\n",
        "      try:\n",
        "        json_exists = os.path.exists(json_name)\n",
        "        if json_exists:\n",
        "          print(\"    > Obs file already exists\")\n",
        "        else:\n",
        "          print(\"    > Grabbing Obs\")\n",
        "          if os.path.exists(\"obs\"):\n",
        "            pass\n",
        "          else:\n",
        "            os.system('mkdir obs')\n",
        "          urlretrieve(obs_url, json_name)\n",
        "        with open(json_name) as json_file:\n",
        "            obs_json = json.load(json_file)\n",
        "            obs_lats = []\n",
        "            obs_lons = []\n",
        "            obs_value = []\n",
        "            obs_elev = []\n",
        "            obs_stid = []\n",
        "            obs_name = []\n",
        "            for stn in obs_json[\"STATION\"]:\n",
        "              # print(stn.encode('utf-8'))\n",
        "              if stn[\"STID\"] is None:\n",
        "                stid = \"N0N3\"\n",
        "              else:\n",
        "                stid = stn[\"STID\"]\n",
        "              #print(f'Processing {region} station {stid}')\n",
        "              name = stn[\"NAME\"]\n",
        "              if stn[\"ELEVATION\"]:\n",
        "                elev = stn[\"ELEVATION\"]\n",
        "              elif stn[\"ELEV_DEM\"]:\n",
        "                elev = stn[\"ELEV_DEM\"]\n",
        "              else:\n",
        "                elev = -999\n",
        "              lat = stn[\"LATITUDE\"]\n",
        "              lon = stn[\"LONGITUDE\"]\n",
        "              stat= None\n",
        "              if \"air_temp_set_1\" in stn['STATISTICS'] and stn['STATISTICS']['air_temp_set_1']:\n",
        "                  if ob_stat in stn['STATISTICS']['air_temp_set_1'] and float(stn[\"LATITUDE\"]) != 0. and float(stn[\"LONGITUDE\"]) != 0.:\n",
        "                      stat = stn['STATISTICS']['air_temp_set_1'][ob_stat]\n",
        "                      obs_stid.append(str(stid))\n",
        "                      obs_name.append(str(name))\n",
        "                      obs_elev.append(int(elev))\n",
        "                      obs_lats.append(float(lat))\n",
        "                      obs_lons.append(float(lon))\n",
        "                      obs_value.append(float(stat))\n",
        "            obs_df = pd.DataFrame()\n",
        "            obs_df[\"stid\"] = obs_stid\n",
        "            obs_df[\"lat\"] = obs_lats\n",
        "            obs_df[\"lon\"] = obs_lons\n",
        "            obs_df[\"maxt\"] = obs_value\n",
        "            obs_csv_name  = \"MaxT_\"+valid_date_end.strftime('%Y%m%d')+\"_\"+dom+\".csv\"\n",
        "            obs_df.to_csv(obs_csv_name)\n",
        "            if element == \"HotMaxT\":\n",
        "              obs_hit = obs_df[(obs_df.maxt > int(threshold_temp))]\n",
        "              obs_miss = obs_df[(obs_df.maxt <= int(threshold_temp))]\n",
        "            elif element == \"ColdMaxT\":\n",
        "              obs_hit = obs_df[(obs_df.maxt < int(threshold_temp))]\n",
        "              obs_miss = obs_df[(obs_df.maxt >= int(threshold_temp))]\n",
        "\n",
        "            obs_hit_lats= obs_hit['lat'].to_numpy()\n",
        "            obs_hit_lons= obs_hit['lon'].to_numpy()\n",
        "            obs_hit_mag= obs_hit['maxt'].to_numpy()\n",
        "            obs_miss_lats= obs_miss['lat'].to_numpy()\n",
        "            obs_miss_lons= obs_miss['lon'].to_numpy()\n",
        "            obs_miss_mag= obs_miss['maxt'].to_numpy()\n",
        "        #plt.scatter(obs_hit_lons, obs_hit_lats, marker=f\"${str(int(round(obs_hit_mag)))}$\", color=\"#4f0615\", edgecolor='none', s=25,transform=latloncrs)\n",
        "        #plt.scatter(obs_hit_lons, obs_hit_lats, marker=get_marker(icon_dict[element][\"Ob_hit\"][\"icon\"]), color=icon_dict[element][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.1, s=35,transform=latloncrs, zorder=3)\n",
        "        #plt.scatter(obs_miss_lons, obs_miss_lats, marker=get_marker(icon_dict[element][\"Ob_miss\"][\"icon\"]), color=icon_dict[element][\"Ob_miss\"][\"color\"], edgecolor='none', s=15,transform=latloncrs, alpha=0.5, zorder=2)\n",
        "        plt.scatter(obs_hit_lons, obs_hit_lats, marker=\"o\", color=icon_dict[element][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.5, s=10, transform=latloncrs, zorder=3)\n",
        "        if plot_ob_misses:\n",
        "          plt.scatter(obs_miss_lons, obs_miss_lats, marker=\"o\", color=\"#000000\", edgecolor=icon_dict[element][\"Ob_miss\"][\"color\"], linewidth=0.5, s=10, transform=latloncrs, alpha=0.3, zorder=2)\n",
        "        print(\"    ✅ So plotted\")\n",
        "      except:\n",
        "        print(\"    ❌ No Obs for you\")\n",
        "\n",
        "    if element == \"HotMaxT\":\n",
        "      title = 'Chance of High Hotter Than '+str(threshold_temp)+\"$^\\circ$F\"\n",
        "    elif element == \"ColdMaxT\":\n",
        "      title = 'Chance of High Colder Than '+str(threshold_temp)+\"$^\\circ$F\"\n",
        "    prodid = 'probmaxt_'+str(valid_date)+'_'+str(threshold_temp)\n",
        "    drawmap(probmaxt_plot, title, prodid, units, CLEVS)\n",
        "\n",
        "def plot_mint_prob():\n",
        "    print(\"Making prob MinT plot...\")\n",
        "    blankmap()\n",
        "    print(\" > Contouring NBM\")\n",
        "    #probmint = nbm\n",
        "    units = r\"%\"\n",
        "    CLEVS = [5,10,20,30,40,50,60,70,80,90,100]\n",
        "    GyBuWh = [\"#4D4D4D\", \"#4D566C\", \"#4E5E84\", \"#4E669A\", \"#4E6DB4\",\n",
        "            \"#5D7CC4\", \"#7D99D0\", \"#9BB2DB\", \"#BCCCE7\", \"#E0EAF5\"]\n",
        "    GyRdYl = [\"#4D4D4D\", \"#724747\", \"#914242\", \"#ad3e3e\", \"#c93939\",\n",
        "            \"#c94b46\", \"#e26e5f\", \"#ec9a7f\", \"#f6c59f\", \"#ffeebd\"]\n",
        "\n",
        "    if \"Cold\" in element:\n",
        "      if \"Dark\" in map_theme:\n",
        "        cmap = GyBuWh\n",
        "      else:\n",
        "        cmap_sel = plt.get_cmap('BuPu')\n",
        "        cmap=[cmap_sel(0.0), cmap_sel(0.1), cmap_sel(0.2), cmap_sel(0.3), cmap_sel(0.4), cmap_sel(0.5), cmap_sel(0.6), cmap_sel(0.7), cmap_sel(0.8), cmap_sel(0.9), cmap_sel(1.0)]\n",
        "    elif \"Hot\" in element:\n",
        "      if \"Dark\" in map_theme:\n",
        "        cmap = GyRdYl\n",
        "      else:\n",
        "        cmap_sel = plt.get_cmap('YlOrRd')\n",
        "        cmap=[cmap_sel(0.0), cmap_sel(0.1), cmap_sel(0.2), cmap_sel(0.3), cmap_sel(0.4), cmap_sel(0.5), cmap_sel(0.6), cmap_sel(0.7), cmap_sel(0.8), cmap_sel(0.9), cmap_sel(1.0)]\n",
        "\n",
        "    probmint_plot = plt.contourf(nbmlons,nbmlats,nbm,CLEVS,colors=cmap, extend='neither',\n",
        "                                 alpha=0.75, transform=latloncrs, transform_first=True, antialiased = True)\n",
        "    #uncomment the following lines if you want contours and labels\n",
        "    #contours = plt.contour(x,y,nbm,CLEVS,colors='silver',linewidths=1., alpha=0.35, transform=latloncrs)\n",
        "    #clabs = plt.clabel(contours, inline=1, fontsize=8, fmt='%1.0f', use_clabeltext=True, colors=font_color)\n",
        "    #plt.setp(clabs, fontweight='bold', path_effects=[PathEffects.withStroke(linewidth=1.4,foreground=face_color)])\n",
        "\n",
        "    if ob_opt:\n",
        "      print(\" > Now looking for obs\")\n",
        "      bbox_str = str(west)+\",\"+str(south)+\",\"+str(east)+\",\"+str(north)\n",
        "      obs_url = \"https://api.synopticlabs.org/v2/stations/statistics?token=\"+synoptic_token+\"&vars=air_temp&start=\"+ \\\n",
        "                valid_date_start.strftime('%Y%m%d')+obs_start_hour+\"&end=\"+valid_date_end.strftime('%Y%m%d')+ \\\n",
        "                obs_end_hour+\"&units=temp%7Cf&within=1440&type=\"+ob_stat+\"&status=active\"+network_string+\"&bbox=\"+bbox_str+\"&fields=latitude,longitude\"\n",
        "      json_name = \"obs/Obs_MinT_\"+valid_date_start.strftime('%Y%m%d')+obs_start_hour+\"_\"+valid_date_end.strftime('%Y%m%d')+obs_end_hour+\"_\"+dom+\".json\"\n",
        "      try:\n",
        "        json_exists = os.path.exists(json_name)\n",
        "        if json_exists:\n",
        "          print(\"    > Obs file already exists\")\n",
        "        else:\n",
        "          print(\"    > Grabbing Obs\")\n",
        "          if os.path.exists(\"obs\"):\n",
        "            pass\n",
        "          else:\n",
        "            os.system('mkdir obs')\n",
        "          urlretrieve(obs_url, json_name)\n",
        "        with open(json_name) as json_file:\n",
        "            obs_json = json.load(json_file)\n",
        "            obs_lats = []\n",
        "            obs_lons = []\n",
        "            obs_value = []\n",
        "            obs_elev = []\n",
        "            obs_stid = []\n",
        "            obs_name = []\n",
        "            for stn in obs_json[\"STATION\"]:\n",
        "              # print(stn.encode('utf-8'))\n",
        "              if stn[\"STID\"] is None:\n",
        "                stid = \"N0N3\"\n",
        "              else:\n",
        "                stid = stn[\"STID\"]\n",
        "              #print(f'Processing {region} station {stid}')\n",
        "              name = stn[\"NAME\"]\n",
        "              if stn[\"ELEVATION\"]:\n",
        "                elev = stn[\"ELEVATION\"]\n",
        "              elif stn[\"ELEV_DEM\"]:\n",
        "                elev = stn[\"ELEV_DEM\"]\n",
        "              else:\n",
        "                elev = -999\n",
        "              lat = stn[\"LATITUDE\"]\n",
        "              lon = stn[\"LONGITUDE\"]\n",
        "              stat= None\n",
        "              if \"air_temp_set_1\" in stn['STATISTICS'] and stn['STATISTICS']['air_temp_set_1']:\n",
        "                  if ob_stat in stn['STATISTICS']['air_temp_set_1'] and float(stn[\"LATITUDE\"]) != 0. and float(stn[\"LONGITUDE\"]) != 0.:\n",
        "                      stat = stn['STATISTICS']['air_temp_set_1'][ob_stat]\n",
        "                      obs_stid.append(str(stid))\n",
        "                      obs_name.append(str(name))\n",
        "                      obs_elev.append(int(elev))\n",
        "                      obs_lats.append(float(lat))\n",
        "                      obs_lons.append(float(lon))\n",
        "                      obs_value.append(float(stat))\n",
        "            obs_df = pd.DataFrame()\n",
        "            obs_df[\"stid\"] = obs_stid\n",
        "            obs_df[\"lat\"] = obs_lats\n",
        "            obs_df[\"lon\"] = obs_lons\n",
        "            obs_df[\"mint\"] = obs_value\n",
        "            obs_csv_name  = \"MinT_\"+valid_date_end.strftime('%Y%m%d')+\"_\"+dom+\".csv\"\n",
        "            obs_df.to_csv(obs_csv_name)\n",
        "            if element == \"ColdMinT\":\n",
        "              obs_hit = obs_df[(obs_df.mint < int(threshold_temp))]\n",
        "              obs_miss = obs_df[(obs_df.mint >= int(threshold_temp))]\n",
        "            elif element == \"HotMinT\":\n",
        "              obs_hit = obs_df[(obs_df.mint > int(threshold_temp))]\n",
        "              obs_miss = obs_df[(obs_df.mint <= int(threshold_temp))]\n",
        "\n",
        "            obs_hit_lats= obs_hit['lat'].to_numpy()\n",
        "            obs_hit_lons= obs_hit['lon'].to_numpy()\n",
        "            obs_hit_mag= obs_hit['mint'].to_numpy()\n",
        "            obs_miss_lats= obs_miss['lat'].to_numpy()\n",
        "            obs_miss_lons= obs_miss['lon'].to_numpy()\n",
        "            obs_miss_mag= obs_miss['mint'].to_numpy()\n",
        "        #plt.scatter(obs_hit_lons, obs_hit_lats, marker=f\"${str(int(round(obs_hit_mag)))}$\", color=\"#4f0615\", edgecolor='none', s=25,transform=latloncrs)\n",
        "        #plt.scatter(obs_hit_lons, obs_hit_lats, marker=get_marker(icon_dict[element][\"Ob_hit\"][\"icon\"]), color=icon_dict[element][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.1, s=35, transform=latloncrs, zorder=3)\n",
        "        #plt.scatter(obs_miss_lons, obs_miss_lats, marker=get_marker(icon_dict[element][\"Ob_miss\"][\"icon\"]), color=icon_dict[element][\"Ob_miss\"][\"color\"], edgecolor='none', s=15, transform=latloncrs, alpha=0.5, zorder=2)\n",
        "        plt.scatter(obs_hit_lons, obs_hit_lats, marker=\"o\", color=icon_dict[element][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.5, s=10, transform=latloncrs, zorder=3)\n",
        "        if plot_ob_misses:\n",
        "          plt.scatter(obs_miss_lons, obs_miss_lats, marker=\"o\", color=\"#000000\", edgecolor=icon_dict[element][\"Ob_miss\"][\"color\"], linewidth=0.5, s=10, transform=latloncrs, alpha=0.3, zorder=2)\n",
        "        print(\"    ✅ So plotted\")\n",
        "      except:\n",
        "        print(\"    ❌ No obs for you\")\n",
        "\n",
        "    if element == \"ColdMinT\":\n",
        "      title = 'Chance of Low Colder Than '+str(threshold_temp)+\"$^\\circ$F\"\n",
        "    elif element == \"HotMinT\":\n",
        "      title = 'Chance of Low Hotter Than '+str(threshold_temp)+\"$^\\circ$F\"\n",
        "\n",
        "    prodid = 'probmint_'+str(valid_date)+'_'+str(threshold_temp)\n",
        "    drawmap(probmint_plot, title, prodid, units, CLEVS)\n",
        "\n",
        "###### PMxWind plots ######################################################################################################\n",
        "\n",
        "def plot_pwind_percentile():\n",
        "  print(\"Making PMxWind percentile plot\")\n",
        "  blankmap()\n",
        "  wind = nbm\n",
        "  WIND_COLS = [\"#103f78\", \"#225ea8\", \"#1d92c0\",\n",
        "                \"#41b7c4\", \"#7fcdbb\", \"#b4d79e\",\n",
        "                \"#dfff9e\", \"#ffffa6\", \"#ffe873\",\n",
        "                \"#ffc400\", \"#ffaa00\", \"#ff5900\",\n",
        "                \"#ff0000\", \"#a80000\", \"#6e0000\",\n",
        "                \"#ffbee8\", \"#ff73de\",]\n",
        "  WIND_LEVS = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 100, 120, 140]\n",
        "\n",
        "  print(\" > Contouring NBM\")\n",
        "  PWIND = plt.contourf(nbmlons,nbmlats,wind,WIND_LEVS,colors=WIND_COLS, extend='neither', alpha=0.75, transform=latloncrs,antialiased = True, transform_first=True)\n",
        "  units = \" mph\"\n",
        "  title = f\"{threshold_percentile}{perc_suffix_dict[threshold_percentile]} Percentile Daily Max Wind {wind_type.capitalize()}\"\n",
        "  prodid = f'pwind{forecast_length}_{threshold_percentile}'\n",
        "\n",
        "  drawmap(PWIND, title, prodid, units, WIND_LEVS)\n",
        "\n",
        "\n",
        "def plot_pwind_prob():\n",
        "  print(\"Making PMxWind probability plot\")\n",
        "  blankmap()\n",
        "  wind = nbm\n",
        "  YlOrBr = [\"#FCF2D5\", \"#F4D1AA\", \"#ECAA79\", \"#E78A51\", \"#E26B32\",\n",
        "            \"#D35625\", \"#BB4A1F\", \"#A13F19\", \"#873213\", \"#64200B\"]\n",
        "\n",
        "  GyOrYl = [\"#5E5E5E\", \"#806551\", \"#986B4A\", \"#AE7042\", \"#C6753D\",\n",
        "            \"#D68342\", \"#DE9750\", \"#E5AB5E\", \"#F7C260\", \"#F9DD83\"]\n",
        "\n",
        "  bilboa = [\"#E5E3E2\", \"#CDC9C0\", \"#BDB18E\", \"#B3956B\", \"#AD765D\",\n",
        "            \"#A65552\", \"#963740\", \"#831F2C\", \"#6B0118\", \"#520000\"]\n",
        "  if \"Dark\" in map_theme:\n",
        "    cmap = GyOrYl\n",
        "  else:\n",
        "    cmap = YlOrBr\n",
        "  CLEVS = [5,10,20,30,40,50,60,70,80,90,100]\n",
        "\n",
        "  print(\" > Contouring NBM\")\n",
        "  PWIND = plt.contourf(nbmlons,nbmlats,wind,CLEVS,colors=cmap, extend='neither', alpha=0.75, transform=latloncrs,antialiased = True, transform_first=True)\n",
        "  units = r\"%\"\n",
        "  title = f\"Chance of Daily Max Wind {wind_type.capitalize()} > {wind_threshold} mph\"\n",
        "  prodid = f'pwind{forecast_length}_gt{wind_threshold}'\n",
        "\n",
        "  if ob_opt:\n",
        "      print(\" > Now looking for obs\")\n",
        "      bbox_str = str(west)+\",\"+str(south)+\",\"+str(east)+\",\"+str(north)\n",
        "      if wind_type == \"gust\":\n",
        "        ob_var = \"wind_gust\"\n",
        "      elif wind_type == \"speed\":\n",
        "        ob_var = \"wind_speed\"\n",
        "      ob_var_set = ob_var + \"_set_1\"\n",
        "      obs_url = \"https://api.synopticlabs.org/v2/stations/statistics?token=\"+synoptic_token+\"&vars=\"+ob_var+\"&start=\"+ \\\n",
        "                valid_date_start.strftime('%Y%m%d')+obs_start_hour+\"&end=\"+valid_date_end.strftime('%Y%m%d')+ \\\n",
        "                obs_end_hour+\"&units=speed|mph&within=1440&type=\"+ob_stat+\"&status=active\"+network_string+\"&bbox=\"+bbox_str+\"&fields=latitude,longitude\"\n",
        "      json_name = \"obs/Obs_\"+ob_var+\"_\"+valid_date_start.strftime('%Y%m%d')+obs_start_hour+\"_\"+valid_date_end.strftime('%Y%m%d')+obs_end_hour+\"_\"+dom+\".json\"\n",
        "      try:\n",
        "        json_exists = os.path.exists(json_name)\n",
        "        if json_exists:\n",
        "          print(\"    > Obs file already exists\")\n",
        "        else:\n",
        "          print(\"    > Grabbing Obs\")\n",
        "          if os.path.exists(\"obs\"):\n",
        "            pass\n",
        "          else:\n",
        "            os.system('mkdir obs')\n",
        "          urlretrieve(obs_url, json_name)\n",
        "        with open(json_name) as json_file:\n",
        "            obs_json = json.load(json_file)\n",
        "            obs_lats = []\n",
        "            obs_lons = []\n",
        "            obs_value = []\n",
        "            obs_elev = []\n",
        "            obs_stid = []\n",
        "            obs_name = []\n",
        "            for stn in obs_json[\"STATION\"]:\n",
        "              # print(stn.encode('utf-8'))\n",
        "              if stn[\"STID\"] is None:\n",
        "                stid = \"N0N3\"\n",
        "              else:\n",
        "                stid = stn[\"STID\"]\n",
        "              #print(f'Processing {region} station {stid}')\n",
        "              name = stn[\"NAME\"]\n",
        "              if stn[\"ELEVATION\"]:\n",
        "                elev = stn[\"ELEVATION\"]\n",
        "              elif stn[\"ELEV_DEM\"]:\n",
        "                elev = stn[\"ELEV_DEM\"]\n",
        "              else:\n",
        "                elev = -999\n",
        "              lat = stn[\"LATITUDE\"]\n",
        "              lon = stn[\"LONGITUDE\"]\n",
        "              stat= None\n",
        "              if ob_var_set in stn['STATISTICS'] and stn['STATISTICS'][ob_var_set]:\n",
        "                  if ob_stat in stn['STATISTICS'][ob_var_set] and float(stn[\"LATITUDE\"]) != 0. and float(stn[\"LONGITUDE\"]) != 0.:\n",
        "                      stat = stn['STATISTICS'][ob_var_set][ob_stat]\n",
        "                      obs_stid.append(str(stid))\n",
        "                      obs_name.append(str(name))\n",
        "                      obs_elev.append(int(elev))\n",
        "                      obs_lats.append(float(lat))\n",
        "                      obs_lons.append(float(lon))\n",
        "                      obs_value.append(float(stat))\n",
        "            obs_df = pd.DataFrame()\n",
        "            obs_df[\"stid\"] = obs_stid\n",
        "            obs_df[\"lat\"] = obs_lats\n",
        "            obs_df[\"lon\"] = obs_lons\n",
        "            obs_df[\"maxwind\"] = obs_value\n",
        "            obs_csv_name  = ob_var+\"_\"+valid_date_end.strftime('%Y%m%d')+\"_\"+dom+\".csv\"\n",
        "            obs_df.to_csv(obs_csv_name)\n",
        "            obs_hit = obs_df[(obs_df.maxwind > float(wind_threshold))]\n",
        "            obs_miss = obs_df[(obs_df.maxwind <= float(wind_threshold))]\n",
        "\n",
        "            obs_hit_lats= obs_hit['lat'].to_numpy()\n",
        "            obs_hit_lons= obs_hit['lon'].to_numpy()\n",
        "            obs_hit_mag= obs_hit['maxwind'].to_numpy()\n",
        "            obs_miss_lats= obs_miss['lat'].to_numpy()\n",
        "            obs_miss_lons= obs_miss['lon'].to_numpy()\n",
        "            obs_miss_mag= obs_miss['maxwind'].to_numpy()\n",
        "\n",
        "        plt.scatter(obs_hit_lons, obs_hit_lats, marker=\"o\", color=icon_dict[nbm_var][\"Ob_hit\"][\"color\"], edgecolor='w', linewidth=0.5, s=10, transform=latloncrs, zorder=3)\n",
        "        if plot_ob_misses:\n",
        "          plt.scatter(obs_miss_lons, obs_miss_lats, marker=\"o\", color=\"#000000\", edgecolor=icon_dict[nbm_var][\"Ob_miss\"][\"color\"], linewidth=0.5, s=10, transform=latloncrs, alpha=0.3, zorder=2)\n",
        "        print(\"    ✅ So plotted\")\n",
        "      except:\n",
        "        print(\"    ❌ No obs for you\")\n",
        "\n",
        "\n",
        "  if lsr_opt:\n",
        "    print(\" > Getting LSRs\")\n",
        "    lsr_start = valid_date_start.strftime('%Y-%m-%dT')\n",
        "    lsr_end = valid_date_end.strftime('%Y-%m-%dT')\n",
        "    lsr_url = \"https://mesonet.agron.iastate.edu/cgi-bin/request/gis/lsr.py?wfo[]=ALL&sts=\"+lsr_start+\"0600Z&ets=\"+lsr_end+\"0600Z&fmt=csv\"\n",
        "    csv_name = \"lsr_\"+lsr_start+\"_\"+lsr_end+\".csv\"\n",
        "\n",
        "    try:\n",
        "      file_exists = os.path.exists(csv_name)\n",
        "      if file_exists:\n",
        "        pass\n",
        "      else:\n",
        "        urlretrieve(lsr_url, csv_name)\n",
        "      LSRs = pd.read_csv(csv_name, usecols= ['LAT','LON','MAG','TYPECODE','TYPETEXT'], on_bad_lines='skip')\n",
        "      LSRs['MAG'] = LSRs.MAG.replace('None','0').astype(float)\n",
        "      LSRs['LAT'] = LSRs.LAT.astype(float)\n",
        "      LSRs['LON'] = LSRs.LON.astype(float)\n",
        "      if wind_type == \"gust\":\n",
        "        lsr_type = 'NON-TSTM WND GST'\n",
        "      elif wind_type == \"speed\":\n",
        "        lsr_type = 'HIGH SUST WINDS'\n",
        "      LSRs_filtered = LSRs[(LSRs.TYPETEXT == lsr_type) & (LSRs.MAG >= float(wind_threshold))]\n",
        "      lsr_lats= LSRs_filtered['LAT'].to_numpy()\n",
        "      lsr_lons= LSRs_filtered['LON'].to_numpy()\n",
        "      lsr_mag= LSRs_filtered['MAG'].to_numpy()\n",
        "      plt.scatter(lsr_lons, lsr_lats, marker=get_marker(icon_dict[nbm_var][\"LSR\"][\"icon\"]), color=icon_dict[nbm_var][\"LSR\"][\"color\"], edgecolor=None, linewidth=0.0, s=15,transform=latloncrs, zorder=3)\n",
        "      print(\"   ✅ ...and plotted!\")\n",
        "    except:\n",
        "      (\"   ❌ No LSRs for you\")\n",
        "\n",
        "  drawmap(PWIND, title, prodid, units, CLEVS)\n",
        "\n",
        "\n",
        "# Running this will make the pretty plot. It might take a few (up to 30) seconds, and will probably spit out some warnings - because it was me who made it after all.\n",
        "matplotlib.rcParams['figure.dpi'] = 200 # make high quality figure - this never worked in the imports\n",
        "if nbm_var == \"PMxMnT\":\n",
        "  get_nbm()\n",
        "  if \"MaxT\" in element:\n",
        "    plot_maxt_prob()\n",
        "  elif \"MinT\" in element:\n",
        "    plot_mint_prob()\n",
        "elif nbm_var == \"PMxWind\":\n",
        "  get_nbm()\n",
        "  if wind_var == \"percentile\":\n",
        "    plot_pwind_percentile()\n",
        "  elif wind_var == \"prob\":\n",
        "    plot_pwind_prob()\n",
        "\n",
        "elif nbm_var == \"PQPF\":\n",
        "  if overview_opt:\n",
        "    from google.colab import widgets\n",
        "    from IPython.utils import io\n",
        "    if ptype == \"qpf\":\n",
        "      if var == \"percentile\":\n",
        "        overview_element_list = [1,10,25,50,75,90,99]\n",
        "        tb = widgets.TabBar([str(i)+perc_suffix_dict[str(i)]+\" Percentile\" for i in overview_element_list])\n",
        "        for p in range(len(overview_element_list)):\n",
        "          threshold_percentile = str(overview_element_list[p])\n",
        "          with io.capture_output() as captured:\n",
        "            get_nbm()\n",
        "          with tb.output_to(p, select=(p < 1)):\n",
        "            plot_qpf_percentile_overview()\n",
        "      if var == \"prob\":\n",
        "        print(\"   ! > Overview option only for PQPF percentiles....for now!\")\n",
        "        '''\n",
        "        overview_element_list = [\"0.01\",\"0.10\",\"0.25\",\"0.50\",\"1.00\",\"2.00\",\"3.00\"]\n",
        "        tb = widgets.TabBar([\"Prob \"+str(i) for i in overview_element_list])\n",
        "        for p in range(len(overview_element_list)):\n",
        "          threshold_percentile = str(overview_element_list[p])\n",
        "          with io.capture_output() as captured:\n",
        "            get_nbm()\n",
        "          with tb.output_to(p, select=(p < 1)):\n",
        "            plot_qpf_percentile_overview()\n",
        "        '''\n",
        "  else:\n",
        "    get_nbm()\n",
        "    if ptype == \"snow\":\n",
        "      if var == \"prob\":\n",
        "        plot_snow_prob()\n",
        "      elif var == \"percentile\":\n",
        "        plot_snow_percentile()\n",
        "      else:\n",
        "        print(\"I'm not sure what you want to do, check your var selection\")\n",
        "    elif ptype == \"qpf\":\n",
        "      if var == \"prob\":\n",
        "        plot_qpf_prob()\n",
        "      elif var == \"percentile\":\n",
        "        plot_qpf_percentile()\n",
        "      else:\n",
        "        print(\"I'm not sure what you want to do, check your var selection\")\n",
        "    elif ptype == \"ice\":\n",
        "      if var == \"prob\":\n",
        "        plot_ice_prob()\n",
        "      elif var == \"percentile\":\n",
        "        plot_ice_percentile()\n",
        "      else:\n",
        "        print(\"I'm not sure what you want to do, check your var selection\")\n",
        "else:\n",
        "  print(\"I'm not sure what you want to do, check your var selection\")"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "Z5xjU68uuyqH",
        "27yZ0aNmu1Up"
      ],
      "provenance": []
    },
    "kernelspec": {
      "display_name": "base",
      "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.9.16"
    },
    "vscode": {
      "interpreter": {
        "hash": "d9efab89797b4f7e4129f7fe7c375038c6a3f1b6c83da7efdea02c4da588d5be"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}