Pirate_HRRR_SM_Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Pirate Weather Example script for SMSL  \n",
    "\n",
    "import os\n",
    "import json\n",
    "import boto3\n",
    "from botocore import UNSIGNED\n",
    "from botocore.client import Config\n",
    "import os\n",
    "import shutil\n",
    "import sys\n",
    "import uuid\n",
    "from urllib.parse import unquote_plus\n",
    "\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "\n",
    "# Plotting tools\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Date processing\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "#os.system('apt install gfortran')\n",
    "\n",
    "# Wgrib2 setup\n",
    "import pywgrib2_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% S3 Client Setup\n"
    }
   },
   "outputs": [],
   "source": [
    "# S3 Client Setup\n",
    "s3_client = boto3.client('s3', config=Config(signature_version=UNSIGNED))\n",
    "s3 = boto3.resource('s3')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Set variables\n"
    }
   },
   "outputs": [],
   "source": [
    "#%% Set download variables\n",
    "bucket          = 'noaa-hrrr-bdp-pds'\n",
    "download_path   = '/tmp'\n",
    "time_in         = '2022-09-16T00:00:00'\n",
    "product         = 'hrrr'\n",
    "\n",
    "# Datetime Setup\n",
    "datetime_RUN        = datetime.strptime(time_in,\"%Y-%m-%dT%H:%M:%S\")\n",
    "fDate          = datetime_RUN.strftime(\"%Y%m%d\")\n",
    "runTime        = \"t\" + datetime_RUN.strftime(\"%H\") + \"z\"\n",
    "\n",
    "# Path Setup\n",
    "download_path_CK    = download_path + '/' + product + '/' + fDate + '/' + runTime + '/' + 'out_' + product + '_chunked.nc'\n",
    "download_path_NC_A  = download_path + '/' + product + '_tmp.nc'\n",
    "\n",
    "if not os.path.exists(download_path + '/' + product):\n",
    "    os.makedirs(download_path + '/' + product)\n",
    "if not os.path.exists(download_path + '/' + product + '/' + fDate):\n",
    "    os.makedirs(download_path + '/' + product + '/' + fDate)\n",
    "if not os.path.exists(download_path + '/' + product + '/' + fDate + '/' + runTime):\n",
    "    os.makedirs(download_path + '/' + product + '/' + fDate + '/' + runTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Setup grid transformation\n"
    }
   },
   "outputs": [],
   "source": [
    "# Setup grid transformation parameters \n",
    "HRRR_grid1 = 'lambert:262.500000:38.500000:38.500000:38.500000'\n",
    "HRRR_grid2 = '237.280472:1799:3000.000000'\n",
    "HRRR_grid3 = '21.138123:1059:3000.000000'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Download Grib File from AWS Open Data\n"
    }
   },
   "outputs": [],
   "source": [
    "# Download Grib File from AWS Open Data\n",
    "# Setup file names\n",
    "grbType   = 'wrfsfc'\n",
    "\n",
    "# Set the range of forecasts to download\n",
    "ncFileRange_A = range(1, 7)\n",
    "\n",
    "for ncFileName in ncFileRange_A:\n",
    "    download_file_pathA  = download_path + '/hrrrh.f' + str(ncFileName).zfill(3) + '.' + grbType + '.grb'\n",
    "    download_path_GB_A   = download_path + '/hrrrh.f' + str(ncFileName).zfill(3) + '.' + grbType + '.grb.earth'\n",
    "\n",
    "    s3_filename = 'hrrr.' + fDate + '/conus/hrrr.' + runTime + '.wrfsfcf' + str(ncFileName).zfill(2) + '.grib2'\n",
    "\n",
    "    # Download the grib file\n",
    "    s3_client.download_file(bucket, s3_filename, download_file_pathA)\n",
    "\n",
    "    #%% Process Downloaded File\n",
    "\n",
    "    # Set variables to read\n",
    "    matchString = (\":(TMP:2 m above ground|CRAIN:surface|CSNOW:surface|\"\n",
    "                   \"CFRZR:surface|PRATE:surface|PRES:surface|CICEP:surface|\"\n",
    "                   \"UGRD:10 m above ground:.*hour fcst|\"\n",
    "                   \"VGRD:10 m above ground:.*hour fcst|\"\n",
    "                   \"VIS:surface|DPT:2 m above ground|TCDC:entire atmosphere|GUST:surface|RH:2 m above ground):\")\n",
    "\n",
    "    # Extract variables and convert winds from grid relative to earth relative\n",
    "    pywgrib2_s.wgrib2([download_file_pathA, '-new_grid_winds', 'earth', '-new_grid_interpolation', 'neighbor', '-match', matchString, '-new_grid', HRRR_grid1, HRRR_grid2, HRRR_grid3, download_path_GB_A])\n",
    "    # Add precipitation\n",
    "    pywgrib2_s.wgrib2([download_file_pathA, '-rewind_init', download_file_pathA, '-new_grid_winds', 'earth', '-new_grid_interpolation', 'neighbor', '-match', 'APCP', '-append','-new_grid', HRRR_grid1, HRRR_grid2, HRRR_grid3, download_path_GB_A, '-quit'])\n",
    "    pywgrib2_s.close(download_path_GB_A)\n",
    "\n",
    "    # Add to NetCDF\n",
    "    pywgrib2_s.wgrib2([download_path_GB_A, '-append', '-netcdf', download_path_NC_A])\n",
    "\n",
    "    pywgrib2_s.close(download_file_pathA)\n",
    "    pywgrib2_s.close(download_path_GB_A)\n",
    "\n",
    "    os.remove(download_file_pathA)\n",
    "    os.remove(download_path_GB_A)\n",
    "\n",
    "    print(ncFileName)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Remove Chunk file if exists\n"
    }
   },
   "outputs": [],
   "source": [
    "# Chunk NetCDF\n",
    "# Remove Chunk file if exists\n",
    "if os.path.isfile(download_path_CK):\n",
    "  os.remove(download_path_CK)\n",
    "if os.path.isfile(download_path  + '/' + product + '/' + fDate + '/' + runTime + '/' + product + '.done'):\n",
    "        os.remove(download_path  + '/' + product + '/' + fDate + '/' + runTime + '/' + product + '.done')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Setup New Datasets\n"
    }
   },
   "outputs": [],
   "source": [
    "# Setup New Datasets\n",
    "chkA    = Dataset(download_path_CK, \"w\")\n",
    "srcA    = Dataset(download_path_NC_A, 'r', format=\"NETCDF3_CLASSIC\")\n",
    "\n",
    "# Copy global attributes all at once via dictionary\n",
    "src = srcA\n",
    "chk = chkA\n",
    "chk.setncatts(srcA.__dict__)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Copy dimensions for srcA and srcB\n"
    }
   },
   "outputs": [],
   "source": [
    "#%% Copy dimensions for source NetCDF\n",
    "for name, dimension in src.dimensions.items():\n",
    "    chk.createDimension(name, (len(dimension) if not dimension.isunlimited() else None))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Copy all file data\n"
    }
   },
   "outputs": [],
   "source": [
    "#%% Copy all file data\n",
    "# This step puts the NetCDF file into time oriented chuncks, as well as adding compression\n",
    "for name, variable in src.variables.items():\n",
    "    if len(variable.dimensions)==3:\n",
    "        if 'PRATE_surface' in name:\n",
    "            x = chk.createVariable(name, variable.datatype, variable.dimensions, chunksizes=[30, 10, 10], zlib=True, least_significant_digit=4, complevel=1)\n",
    "        else:\n",
    "            x = chk.createVariable(name, variable.datatype, variable.dimensions, chunksizes=[30, 10, 10], zlib=True, least_significant_digit=1, complevel=1)\n",
    "    else:\n",
    "        x = chk.createVariable(name, variable.datatype, variable.dimensions)\n",
    "\n",
    "    chk[name][:] = src[name][:]\n",
    "    # copy variable attributes all at once via dictionary\n",
    "    for ncattr in src[name].ncattrs():\n",
    "        if ncattr != '_FillValue':\n",
    "            chk[name].setncattr(ncattr, src[name].getncattr(ncattr))\n",
    "\n",
    "src.close()\n",
    "chk.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Extract data for a given lat lon from NetCDF\n"
    }
   },
   "outputs": [],
   "source": [
    "ChunkNCfile = Dataset(download_path_CK, 'r')\n",
    "\n",
    "# Read Latitude, Longitude, and Times points in file\n",
    "lats_hrrr = ChunkNCfile['latitude'][:, :]\n",
    "lons_hrrr = ChunkNCfile['longitude'][:, :]\n",
    "times_hrrr = ChunkNCfile['time'][:]\n",
    "\n",
    "# Convert timestamps to datetimes\n",
    "datetime_hrrr = []\n",
    "for i in range(0, len(times_hrrr)):\n",
    "    datetime_hrrr.append(datetime.fromtimestamp(times_hrrr[i]))\n",
    "\n",
    "# Set desired lat long\n",
    "##### NOTE #####\n",
    "# Must be in the continental United States or southern Canada\n",
    "Point_Lat = 45.4215\n",
    "Point_Lon = -75.6972+360\n",
    "\n",
    "# Find cloest point\n",
    "abslat = np.abs(lats_hrrr-Point_Lat)\n",
    "abslon = np.abs(lons_hrrr-Point_Lon)\n",
    "c = np.maximum(abslon, abslat)\n",
    "x_hrrr, y_hrrr = np.where(c == np.min(c))\n",
    "\n",
    "# Exract and print temperature forecast!\n",
    "vard = ChunkNCfile['TMP_2maboveground'][:, x_hrrr, y_hrrr]\n",
    "\n",
    "# Convert to degrees Celsius\n",
    "print(vard-272.15)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%% Plot the extracted data\n"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "plt.plot(datetime_hrrr, vard[:, 0, 0]-272.15)\n",
    "plt.xlabel('Datetime')\n",
    "plt.ylabel('Temperature [degC]')\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "instance_type": "ml.m5.large",
  "kernelspec": {
   "display_name": "default:Python",
   "language": "python",
   "name": "conda-env-default-py"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}