Plot_NEXRAD_mosaics.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/tjturnage/radar/blob/main/Plot_NEXRAD_mosaics.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BxsBvynDTsJy"
      },
      "source": [
        "<font size=\"+4\" color=\"green\"><b>Plot NEXRAD Reflectivity Mosaics</b></font>  \n",
        "<font size=\"+1\" color=\"gray\"><i>updated July 30, 2024</i></font>  \n",
        "<hr>  \n",
        "\n",
        "This allows the user to plot nexrad mosaic reflectivity images described at:\n",
        "https://mesonet.agron.iastate.edu/docs/nexrad_mosaic/\n",
        "<br />\n",
        "<br />\n",
        "The user can define the following:\n",
        "- whether to use high or low resolution plots (low-res goes back to 1995)\n",
        "- form-based selection of time range to plot\n",
        "- time zone to include in the label\n",
        "- plotting interval (that is, make a plot for every 5, 10, 15, etc., minutes)\n",
        "- lat/lon box containing the plot\n",
        "- one of two available reflectivity colormaps (additional maps available on request)\n",
        "- which cities to plot (if any) and where to place the labels\n",
        "- whether to plot major roads\n",
        "- whether to plot county outlines\n",
        "- figure size and plot aspect ratio (usually don't need to change)\n",
        "- whether to create an animated gif\n",
        "<br/> <br/>\n",
        "\n",
        "A file called \"radar_images.zip\" will be created. This file can be downloaded and then unzipped on your local machine to extract the image files.\n",
        "\n",
        "\n",
        "---\n",
        "Questions or bugs? Please contact me at thomas.turnage@noaa.gov\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iR0dkaxB0U2q",
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c2616459-883f-4c9b-97b6-d3ddea3750a4"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "2024-07-15 21:00:00 2024-07-16 09:00:00\n"
          ]
        }
      ],
      "source": [
        "# @title <font size=\"+3\" color=\"green\">Select your configuration</font> ...\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Select the type of plot</font></b><br>Recommendation: Choose <b>\"low res\"</b> for multi-state plots</i></b></font>\n",
        "from datetime import datetime, timedelta\n",
        "import pytz\n",
        "\n",
        "\n",
        "plot_type = \"low res (available since 1995)\" # @param [\"low res (available since 1995)\", \"high res (available since Nov2010)\"]\n",
        "\n",
        "file_type = 'n0r'\n",
        "if plot_type == \"high res (available since Nov2010)\":\n",
        "  file_type = 'n0q'\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Enter the start time (in UTC) for the plots</i></b></font>\n",
        "start_date = \"2024-07-15\" # @param {type:\"date\"}\n",
        "start_date = start_date.replace('-', '')\n",
        "start_hour = \"21\" # @param [\"00\",\"01\",\"02\",\"03\",\"04\",\"05\",\"06\",\"07\",\"08\",\"09\",\"10\",\"11\",\"12\",\"13\",\"14\",\"15\",\"16\",\"17,\",\"18\",\"19\",\"20\",\"21\",\"22\",\"23\"]\n",
        "start_minute = \"00\" # @param [\"00\", \"05\", \"10\", \"15\", \"20\", \"25\", \"30\", \"35\", \"40\", \"45\", \"50\", \"55\"]\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Enter the duration in hours and minutes.</b></font><br>Example: 02 hrs, 20 mins = 2:20 duration</font>\n",
        "hours_duration = \"12\" # @param [\"00\",\"01\",\"02\",\"03\",\"04\",\"05\",\"06\",\"07\",\"08\",\"09\",\"10\",\"11\",\"12\"]\n",
        "plus_minutes = \"00\" # @param [\"00\",\"10\",\"20\",\"30\",\"40\",\"50\"]\n",
        "if start_date < '20101113' and file_type == 'n0q':\n",
        "  print(\"High res not available for this time! ...\")\n",
        "  print(\"... selecting low res instead ...\")\n",
        "  file_type = 'n0r'\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Select time zone for labelling your figures</i></b></font>\n",
        "time_zone = \"CDT\" # @param [\"UTC\", \"EST\", \"EDT\", \"CST\", \"CDT\", \"MST\", \"MDT\", \"PST\", \"PDT\"]\n",
        "time_zone_dict = {\"UTC\": 0, \"EST\":5, \"EDT\":4, \"CST\":6, \"CDT\":5, \"MST\":7, \"MDT\":6, \"PST\":8, \"PDT\":7 }\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Define time interval in minutes between plots</i></b></font>\n",
        "time_step = \"10\" # @param [5, 10, 15, 30, 45, 60]\n",
        "time_step = int(time_step)\n",
        "\n",
        "start_datetime_string = f'{start_date} {start_hour}:{start_minute}:00'\n",
        "\n",
        "start_datetime = datetime.strptime(start_datetime_string, '%Y%m%d %H:%M:%S')\n",
        "end_datetime = start_datetime + timedelta(hours=int(hours_duration), minutes=int(plus_minutes))\n",
        "first_plot_time = f'{start_date}_{start_hour}{start_minute}'\n",
        "last_plot_time = datetime.strftime(end_datetime,'%Y-%m-%d_%H%M')\n",
        "#last_plot_time = f'{end_date}_{end_hour}{end_minute}'\n",
        "print(start_datetime, end_datetime)\n",
        "#if first_plot_time >= last_plot_time:\n",
        "#  print(\"... Your start date/time must be before your end date/time !! ...\")\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Enter lat/lon values to define your plotting box</i></b></font>\n",
        "lat_min = 38 # @param {type:\"raw\"}\n",
        "lat_max = 44.0 # @param {type:\"raw\"}\n",
        "lon_min = -94.5 # @param {type:\"raw\"}\n",
        "lon_max = -83.0 # @param {type:\"raw\"}\n",
        "\n",
        "def make_float(val):\n",
        "  try:\n",
        "    return float(val)\n",
        "  except:\n",
        "    print(\"...... WARNING !! ......................\\n\")\n",
        "    print(f\"Can't convert {val} to float!\\n\")\n",
        "    print(\"...... Check values and try again ......\\n\")\n",
        "    return\n",
        "\n",
        "#for coor in (lat_min, lat_max, lon_min, lon_max):\n",
        "\n",
        "lat_min = make_float(lat_min)\n",
        "lat_max = make_float(lat_max)\n",
        "lon_min = make_float(lon_min)\n",
        "lon_max = make_float(lon_max)\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Pick which features to include in the plot</i></b></font>\n",
        "plot_counties = True # @param {type:\"boolean\"}\n",
        "plot_roads = False # @param {type:\"boolean\"}\n",
        "# @markdown <font size=\"+1\" color=\"blue\"><b><i>(Note: if plotting cities, you'll need to create your own dictionary)</i></b></font>\n",
        "plot_cities = False # @param {type:\"boolean\"}\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Check box for a \"snow-friendly\" color table</i></b></font>\n",
        "\n",
        "snow_colormap = False # @param {type:\"boolean\"}\n",
        "if snow_colormap:\n",
        "  cmap_type = \"winter\"\n",
        "else:\n",
        "  cmap_type = \"wdtd\"\n",
        "\n",
        "plot_key = f'{file_type}_{cmap_type}'\n",
        "plot_dict = {'n0r_wdtd':{'cmap':'z_lowres','vmin':-30, 'vmax': 75, 'mult':5, 'subt':30.0},\n",
        "             'n0q_wdtd':{'cmap':'wdtd_z','vmin':-32.5, 'vmax': 90, 'mult':0.5, 'subt':32.5},\n",
        "             'n0r_winter':{'cmap':'winter_z','vmin':-30, 'vmax': 75, 'mult':5, 'subt':30.0},\n",
        "             'n0q_winter':{'cmap':'winter_z','vmin':-32.5, 'vmax': 90, 'mult':0.5, 'subt':32.5}\n",
        "             }\n",
        "\n",
        "plot_values = plot_dict[plot_key]\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>If needed, adjust for undesired map stretching</b></font><br>Edit one or more of the following values - may require some testing</i></b></font>\n",
        "# @markdown <font color=\"blue\"><ol><li>Adjust your lat/lon boundary ranges</li>\n",
        "# @markdown <li>Adjust the height and width dimensions of the figure.</li>\n",
        "# @markdown <li>Adjust the plot aspect ratio ... >1 shrinks horizontal scale and\n",
        "# @markdown vice versa.</li></ol></h4>\n",
        "\n",
        "fig_height = 8 # @param {type:\"raw\"}\n",
        "fig_height = int(fig_height)\n",
        "fig_width = 5 # @param {type:\"raw\"}\n",
        "fig_width = int(fig_width)\n",
        "plot_aspect_ratio = 1.25 # @param {type:\"raw\"}\n",
        "plot_aspect_ratio = float(plot_aspect_ratio)\n",
        "\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>Check to create animated gif</b></font><br>This will take time to load the software package ... creates a file called \"animation.gif\"</i>\n",
        "create_animation = True # @param {type:\"boolean\"}\n",
        "# @markdown <font size=\"+2\" color=\"blue\"><b><i>This sets time between frames in hundredths of a second</i></font></h4>\n",
        "animation_speed =\"20\" # @param [5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100]\n",
        "animation_speed = int(animation_speed)\n",
        "\n",
        "cities = {\n",
        "'Grand Rapids':{'lat':42.963,'lon':-85.668, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Muskegon':{'lat':43.234,'lon':-86.248, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Holland':{'lat':42.7875,'lon':-86.1089, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Lawton':{'lat':42.1675,'lon':-85.8480, 'lat_offset':-0.05, 'lon_offset':0.05},\n",
        "'Marshall':{'lat':42.2722,'lon':-84.9633, 'lat_offset':-0.08, 'lon_offset':0.04},\n",
        "'Battle Creek':{'lat':42.3211,'lon':-85.1800, 'lat_offset':0.04, 'lon_offset':-0.06},\n",
        "'South Haven':{'lat':42.4031,'lon':-86.2736, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Lansing':{'lat':42.7325,'lon':-84.5556, 'lat_offset':0.00, 'lon_offset':0.06},\n",
        "'Kalamazoo':{'lat':42.292,'lon':-85.587, 'lat_offset':0.04, 'lon_offset':-0.14},\n",
        "'Ludington':{'lat':43.955,'lon':-86.4525, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Big Rapids':{'lat':43.698,'lon':-85.4836, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Mount Pleasant':{'lat':43.5978,'lon':-84.7675, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "#'MC109':{'lat':43.0534,'lon':-85.9080, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "#'MC158':{'lat':43.0923,'lon':-86.0977, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "#'MC159':{'lat':43.4101,'lon':-86.3146, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "           }\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pOgLNG5O9bHF",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title <font size=\"+3\" color=\"green\">If you want to customize your own cities dictionary, click \"Show code\" below</font>\n",
        "#   You'll want to test the best offset position for each of the labels.\n",
        "#   Positive \"lat_offset\" shifts the label north by that many degrees\n",
        "#   Positive \"lon_offest\" shifts the label east by that many degrees\n",
        "\n",
        "cities = {\n",
        "'Grand Rapids':{'lat':42.963,'lon':-85.668, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Muskegon':{'lat':43.234,'lon':-86.248, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Holland':{'lat':42.7875,'lon':-86.1089, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'South Haven':{'lat':42.4031,'lon':-86.2736, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Lansing':{'lat':42.7325,'lon':-84.5556, 'lat_offset':0.00, 'lon_offset':0.06},\n",
        "'Kalamazoo':{'lat':42.292,'lon':-85.587, 'lat_offset':0.08, 'lon_offset':0.04},\n",
        "'Portage':{'lat':42.20,'lon':-85.59, 'lat_offset':-0.08, 'lon_offset': 0.04},\n",
        "'Ludington':{'lat':43.955,'lon':-86.4525, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Big Rapids':{'lat':43.698,'lon':-85.4836, 'lat_offset':0.04, 'lon_offset':0.06},\n",
        "'Mount Pleasant':{'lat':43.5978,'lon':-84.7675, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "'Portage':{'lat':42.20,'lon':-85.59, 'lat_offset':-0.08, 'lon_offset': 0.04},\n",
        "#'MC109':{'lat':43.0534,'lon':-85.9080, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "#'MC158':{'lat':43.0923,'lon':-86.0977, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "#'MC159':{'lat':43.4101,'lon':-86.3146, 'lat_offset':0.06, 'lon_offset':0.06},\n",
        "           }\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "aTiQtbh-0P9A"
      },
      "outputs": [],
      "source": [
        "# @title <font size=\"+3\" color=\"green\"> Install the needed GIS software, prepare colormaps and GIS files</font>\n",
        "!pip install rasterio &> /dev/null\n",
        "!pip install cartopy &> /dev/null\n",
        "import sys\n",
        "import os\n",
        "import subprocess\n",
        "import glob\n",
        "from os.path import exists\n",
        "from matplotlib.colors import LinearSegmentedColormap\n",
        "import matplotlib.colors as mcolors\n",
        "import matplotlib as mpl\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import rasterio\n",
        "import urllib.request\n",
        "import requests\n",
        "from bs4 import BeautifulSoup\n",
        "import zipfile\n",
        "import cartopy.crs as ccrs\n",
        "import cartopy.feature as cfeature\n",
        "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
        "import cartopy.io.shapereader as shpreader\n",
        "import xarray as xr\n",
        "from datetime import datetime\n",
        "from datetime import timedelta\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\")\n",
        "\n",
        "os.chdir('/content')\n",
        "\n",
        "def create_and_register_cmap(cmap_name, original_colors, position, rgb_input=True) -> None:\n",
        "    \"\"\"\n",
        "    1) Converts list of color values from 8-bit (0-255 scale)\n",
        "       to 0-1 decimal scale.\n",
        "\n",
        "    2) The created final_color_list is combine with position list\n",
        "       to create a matplotlib colormap.\n",
        "\n",
        "    3) Registers the new colormap using the supplied cmap_name.\n",
        "\n",
        "    Input\n",
        "    ----------\n",
        "             cmap_name : string\n",
        "                         Name of the cmap you want to use and register\n",
        "\n",
        "      original_ colors : list\n",
        "                         8-bit tuples containing RGB values\n",
        "                         example: [(127,256,192), (125,125,125)]\n",
        "\n",
        "              position : list\n",
        "                         node locations for color break points on\n",
        "                         a 0 to 1 scale\n",
        "\n",
        "            rgb_input : boolean\n",
        "                        True: Color conversion from RGB to performed.\n",
        "                       False: Color conversion from RGB is skipped.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    None\n",
        "\n",
        "    \"\"\"\n",
        "    if rgb_input:\n",
        "        final_color_list = []\n",
        "        for color in original_colors:\n",
        "            converted_color_group = []\n",
        "            for original_channel_value in color:\n",
        "                new_channel_value = float(original_channel_value)/255\n",
        "                converted_color_group.append(float(format(new_channel_value, '.3f')))\n",
        "            converted_color_tuple = tuple(converted_color_group)\n",
        "            final_color_list.append(converted_color_tuple)\n",
        "    else:\n",
        "        final_color_list = original_colors\n",
        "\n",
        "    this_cmap = LinearSegmentedColormap.from_list(cmap_name, list(zip(position, final_color_list)))\n",
        "    mpl.colormaps.register(cmap=this_cmap, force=True)\n",
        "\n",
        "\n",
        "###############################################################################\n",
        "\n",
        "winter_rgbs = [(0, 0, 0), (70, 70, 70), (210, 210, 210),\n",
        "(0, 120, 255), (0, 0, 60), (225, 225, 0), (220, 100, 0), (220, 100, 0)]\n",
        "\n",
        "winter_nodes = [0.0, 0.20, 0.38, 0.42, 0.48, 0.52, 0.60, 1.0]\n",
        "\n",
        "create_and_register_cmap('winter_z', winter_rgbs, winter_nodes)\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "\n",
        "wdtd_z_rgbs =[ (0, 0, 0), (0, 0, 0), (0, 0, 0), (130, 130, 130),\n",
        " (95, 189, 207), (57, 201, 105), (57, 201, 105), (0, 40, 0),\n",
        " (9, 94, 9), (255, 207, 0), (255, 207, 0), (255, 133, 0),\n",
        " (255, 0, 0), (89, 0, 0), (255, 245, 255), (225, 11, 227),\n",
        " (164, 0, 247), (99, 0, 214), (5, 221, 224), (58, 103, 181),\n",
        " (255, 255, 255),(255, 255, 255)]\n",
        "\n",
        "def create_nodes(refs,denom,below_zero):\n",
        "  nodes = []\n",
        "  for r in refs:\n",
        "    ref = (r + below_zero)/denom\n",
        "    nodes.append(ref)\n",
        "  return nodes\n",
        "\n",
        "refs = [-32, -25.1, -25,14.99, 15,19.99, 20,34.99, 35,39.99, 40,49.99,\n",
        "        50,59.99, 60,69.99, 70,74.99, 75,79.99, 80,90]\n",
        "\n",
        "wdtd_z_nodes = create_nodes(refs,122,32)\n",
        "\n",
        "refs = [-30, -25.1, -25, 14.99, 15,19.99, 20,34.99, 35,39.99, 40,49.99,\n",
        "        50,59.99, 60,69.99, 70,74.99, 75]\n",
        "\n",
        "low_nodes = create_nodes(refs,105,30)\n",
        "\n",
        "create_and_register_cmap('wdtd_z', wdtd_z_rgbs, wdtd_z_nodes)\n",
        "create_and_register_cmap('z_lowres', wdtd_z_rgbs[:-3], low_nodes)\n",
        "\n",
        "\n",
        "###############################################################################\n",
        "\n",
        "#scan this webpage to determine the most recent NWS County shapefile available\n",
        "url = \"https://www.weather.gov/gis/Counties\"\n",
        "\n",
        "response = requests.get(url)\n",
        "soup = BeautifulSoup(response.content, \"html.parser\")\n",
        "zip_links = soup.find_all(\"a\", href=lambda href: href.endswith(\".zip\"))\n",
        "if zip_links:\n",
        "    zip_url = zip_links[-1][\"href\"]\n",
        "    counties_url = os.path.join('https://www.weather.gov',str(zip_url)[1:])\n",
        "    zip_file_name = counties_url.split(\"/\")[-1]\n",
        "    zip_file_root = zip_file_name[:-4]\n",
        "    shapefile_name = f'{zip_file_root}.shp'\n",
        "else:\n",
        "    print(\"No zip files found on the website.\")\n",
        "\n",
        "if not exists(zip_file_name):\n",
        "  urllib.request.urlretrieve(counties_url, \"counties.zip\")\n",
        "  with zipfile.ZipFile('counties.zip', 'r') as zip_ref:\n",
        "      zip_ref.extractall('.')\n",
        "\n",
        "try:\n",
        "  COUNTY_SHAPEFILE\n",
        "except NameError:\n",
        "  #reader = shpreader.Reader(\"tl_2019_us_county.shp\")\n",
        "  reader = shpreader.Reader(shapefile_name)\n",
        "  features = list(reader.geometries())\n",
        "  COUNTY_SHAPEFILE = cfeature.ShapelyFeature(features, ccrs.PlateCarree())\n",
        "\n",
        "# download and read road shapefile only if plot_roads=True\n",
        "if plot_roads:\n",
        "  road_shapefile_name = \"tl_2019_us_primaryroads.shp\"\n",
        "  road_url = \"https://www2.census.gov/geo/tiger/TIGER2019/PRIMARYROADS/tl_2019_us_primaryroads.zip\"\n",
        "  if not exists(road_shapefile_name):\n",
        "    urllib.request.urlretrieve(road_url, \"roads.zip\")\n",
        "    with zipfile.ZipFile('roads.zip', 'r') as zip_roads:\n",
        "      zip_roads.extractall('.')\n",
        "  try:\n",
        "    ROAD_SHAPEFILE\n",
        "  except NameError:\n",
        "    #reader = shpreader.Reader(\"tl_2019_us_primary_roads.shp\")\n",
        "    road_reader = shpreader.Reader(road_shapefile_name)\n",
        "    features = list(road_reader.geometries())\n",
        "    ROAD_SHAPEFILE = cfeature.ShapelyFeature(features, ccrs.PlateCarree())\n",
        "\n",
        "\n",
        "\n",
        "plt.rcParams['font.size'] = 10\n",
        "plt.rcParams['font.weight'] = 'bold'\n",
        "plt.rcParams['lines.markersize'] = 6\n",
        "plt.rcParams['axes.titlesize'] = 16\n",
        "plt.rcParams['text.color'] = 'white'\n",
        "plt.rcParams['ytick.labelsize'] = 12\n",
        "plt.rcParams['legend.fontsize'] = 10\n",
        "plt.rcParams['scatter.marker'] = 'o'\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0JU4kGjidNj4",
        "cellView": "form",
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "# @title <font size=\"+3\" color=\"green\">Create your images!</font>\n",
        "import matplotlib as mpl\n",
        "\n",
        "def create_valid_times(start_dt, end_dt):\n",
        "\n",
        "  #start_time = datetime.strptime(initial_time, \"%Y%m%d_%H%M\")\n",
        "  #end_time = datetime.strptime(final_time, \"%Y%m%d_%H%M\")\n",
        "\n",
        "  valid_times = []\n",
        "  i = 0\n",
        "  this_time = start_dt + timedelta(minutes=time_step*i)\n",
        "  valid_times.append(this_time)\n",
        "  while this_time < end_dt:\n",
        "    i += 1\n",
        "    this_time = start_dt + timedelta(minutes=time_step*i)\n",
        "    valid_times.append(this_time)\n",
        "\n",
        "  return valid_times\n",
        "\n",
        "valid_times = create_valid_times(start_datetime,end_datetime)\n",
        "\n",
        "\n",
        "# ------------- Download pngs, create plots, save images -------------\n",
        "\n",
        "URL_BASE = \"https://mesonet.agron.iastate.edu/archive/data/\"\n",
        "\n",
        "os.chdir('/content')\n",
        "# create fresh images directory to store new images to\n",
        "result = subprocess.run([\"rm -rf images ; mkdir images\"], shell=True, capture_output=True, text=True)\n",
        "#fout = open('output.txt', 'w')\n",
        "#fout.write(f'dts,station,ref\\n')\n",
        "for plot_time in valid_times:\n",
        "  local_time = plot_time - timedelta(hours=time_zone_dict[time_zone])\n",
        "\n",
        "  file_dirs = datetime.strftime(plot_time,\"%Y/%m/%d/GIS/uscomp/\")\n",
        "  file_timestamp = datetime.strftime(plot_time,\"_%Y%m%d%H%M\")\n",
        "  url_path = f'{URL_BASE}{file_dirs}'\n",
        "  file_root = f'{file_type}{file_timestamp}'\n",
        "  downloaded_png = f'{file_root}.png'\n",
        "  png_url = f'{url_path}{file_root}.png'\n",
        "  wld_url = f'{url_path}{file_root}.wld'\n",
        "  if not exists(downloaded_png):\n",
        "    map = urllib.request.urlretrieve(png_url,downloaded_png)\n",
        "\n",
        "  with rasterio.open(downloaded_png) as src:\n",
        "    # Get the image data and metadata\n",
        "    image_data = src.read()\n",
        "\n",
        "    #slice (x,y) cartesian array out of current (1,x,y) array\n",
        "    image_new = image_data[0,:,:]\n",
        "\n",
        "    # convert 8-bit (0-255) values to reflectivity values as described at\n",
        "    # https://mesonet.agron.iastate.edu/GIS/rasters.php?rid=2\n",
        "\n",
        "    image_new = (image_new * plot_values['mult']) - plot_values['subt']\n",
        "\n",
        "\n",
        "    #flip data vertically to move origin from upper left to lower left\n",
        "    image_flipped = np.flipud(image_new)\n",
        "\n",
        "    # determine images dimensions\n",
        "    xlen = np.shape(image_flipped)[1]\n",
        "    ylen = np.shape(image_flipped)[0]\n",
        "\n",
        "    # create list of lat/lon values based on dimension lengths\n",
        "    x = list(np.linspace(-126.0, -65.0, num=xlen, retstep=True)[0])\n",
        "    y = list(np.linspace(23.0, 50.0, num=ylen, retstep=True)[0])\n",
        "\n",
        "    # Create xarray add latitude/longitude dimensions using lists above\n",
        "    coords = {'longitude': x, 'latitude': y}\n",
        "    data = xr.DataArray(image_flipped, coords=coords, dims=['latitude', 'longitude'])\n",
        "    data.attrs[\"units\"] = \"dBZ\"\n",
        "    data.attrs[\"long_name\"] = \"Reflectivity\"\n",
        "\n",
        "    # Determine region you want to plot\n",
        "    subset = data.sel(latitude=slice(lat_min, lat_max), longitude=slice(lon_min, lon_max))\n",
        "\n",
        "    fig, ax = plt.subplots(1,1,figsize=(fig_height,fig_width),layout='constrained',subplot_kw={'projection': ccrs.PlateCarree()})\n",
        "    subset.plot(ax=ax, cmap=plot_values['cmap'], vmin=plot_values['vmin'], vmax=plot_values['vmax'])\n",
        "\n",
        "\n",
        "    ax.set_aspect(plot_aspect_ratio)\n",
        "    if plot_counties:\n",
        "      ax.add_feature(COUNTY_SHAPEFILE, facecolor='none', edgecolor='#cccccc', linewidth=0.7)\n",
        "\n",
        "    if plot_roads:\n",
        "      ax.add_feature(ROAD_SHAPEFILE, facecolor='none', edgecolor='cyan', linewidth=0.5)\n",
        "    ax.add_feature(cfeature.STATES, facecolor='none', edgecolor='white', linewidth=1.3)\n",
        "    ax.set_title(datetime.strftime(local_time,f'%b %d, %Y -- %I:%M %p {time_zone}'), color='black')\n",
        "    if time_zone == 'UTC':\n",
        "      ax.set_title(datetime.strftime(local_time,f'%b %d, %Y -- %H:%M {time_zone}'), color='black')\n",
        "    #plt.tight_layout()\n",
        "    target_fname = f'images/cropped_{downloaded_png}'\n",
        "    if plot_cities:\n",
        "      for key in cities:\n",
        "        c = cities[key]\n",
        "        lat,lon = c['lat'], c['lon']\n",
        "        if ((lat > lat_min) & (lat < lat_max) & (lon > lon_min) & (lon < lon_max)):\n",
        "          if 'MC' in key:\n",
        "            plt.plot([lon], [lat], color='#FFAAAA', linewidth=1, marker='o', markersize=4, zorder=10)\n",
        "            plt.text(lon + c['lon_offset'], lat + c['lat_offset'], key, color='#FFAAAA', zorder=11)\n",
        "          else:\n",
        "            plt.plot([lon], [lat], color='#FFFF00', linewidth=1, marker='o', markersize=6, zorder=10)\n",
        "            plt.plot([lon], [lat], color='white', linewidth=1, marker='o', markersize=4, zorder=10)\n",
        "            plt.text(lon + c['lon_offset'], lat + c['lat_offset'], key, color='#FFFFFF', zorder=11)\n",
        "    plt.savefig(target_fname,bbox_inches='tight')\n",
        "\n",
        "\n",
        "import glob\n",
        "import os\n",
        "\n",
        "if create_animation:\n",
        "  !apt install imagemagick &> /dev/null\n",
        "  os.chdir('/content/images')\n",
        "  cmd = f'convert -delay 25 -loop 0 *.png /content/animation.gif'\n",
        "  os.system(cmd)\n",
        "\n",
        "# tar images files together so you can download single file\n",
        "os.chdir('/content/')\n",
        "with zipfile.ZipFile('radar_images.zip', 'w') as f:\n",
        "    for file in glob.glob('images/*'):\n",
        "        f.write(file)\n",
        "\n",
        "\n",
        "with zipfile.ZipFile('/content/images/images.zip', 'r') as zip_ref:\n",
        "    zip_ref.extractall('/content/images')\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}