import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from netCDF4 import num2date
import numpy as np
import xarray as xr
from siphon.catalog import TDSCatalog
from datetime import datetime
import datetime as dt
from xarray.backends import NetCDF4DataStore
import cartopy
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.ndimage import gaussian_filter
import metpy.calc as mpcalc
import numpy.ma as ma
from metpy.units import units
import scipy.ndimage as ndimage
from metpy.plots import USCOUNTIES
import matplotlib.patches as mpatches

def mkdir_p(mypath):
    '''Creates a directory. equivalent to using mkdir -p on the command line'''

    from errno import EEXIST
    from os import makedirs,path

    try:
        makedirs(mypath)
    except OSError as exc: # Python >2.5
        if exc.errno == EEXIST and path.isdir(mypath):
            pass
        else: raise

startTime=datetime.now()

year = startTime.year

if startTime.month <10:
    month = '0'+str(startTime.month)
else:
    month = str(startTime.month)

if startTime.day <10:
    day = '0'+str(startTime.day)
else:
    day = str(startTime.day)

if startTime.hour <10:
    hour = '0'+str(startTime.hour)
else:
    hour = str(startTime.hour)

mdate = str(year)+str(month)+str(day)

def get_init_hr(hour):
    if int(hour) <3:
        init_hour = '00'
    elif int(hour) <9:
        init_hour = '06'
    elif int(hour) <15:
        init_hour = '12'
    elif int(hour) <21:
        init_hour = '18'
    else:
        init_hour = '00'
    return(init_hour)
init_hour = get_init_hr(hour)
url = 'http://nomads.ncep.noaa.gov:80/dods/hrrr/hrrr'+mdate+'/hrrr_sfc.t'+get_init_hr(hour)+'z'
#url='http://nomads.ncep.noaa.gov:80/dods/hrrr/hrrr20201231/hrrr_sfc.t00z'
print(url)

# Create new directory
output_dir = str(year)+str(month)+str(day)+'_'+str(init_hour)+'00'
#output_dir = '20201231_0000'
mkdir_p(output_dir)
mkdir_p(output_dir+'/HRRR_test')
#Parse data using MetPy
ds = xr.open_dataset(url)
init_hr = dt.datetime(int(year),int(month),int(day),int(init_hour))
times = ds['tmp2m'].metpy.time
init_time = ds['time'][0]

lats = np.arange(25,55,0.25)
lons = np.arange(260,310,0.25)

for i in range(1,49):
    fc_hr = init_hr+dt.timedelta(hours=1*i)
    forecast_hour = times[0].values

    data = ds.metpy.parse_cf()
    data = data.isel(time=i)

    #Rename variables to useful things
    data = data.rename({
        'cfrzrsfc':'catice',
        'cicepsfc':'catsleet',
        'crainsfc':'catrain',
        'csnowsfc':'catsnow',
        'tcdcclm':'tcc',
        'tmpprs': 'temperature',
        'ugrd10m': 'u',
        'vgrd10m': 'v',
        'mslmamsl':'mslp',
        'tmp2m':'sfc_temp',
        'dpt2m':'sfc_td',
        'refcclm':'radar',
        'apcpsfc':'qpf',
        'capesfc':'cape',
        'gustsfc':'sfcgust',
        'hcdchcll':'high_cloud',
        'mcdcmcll':'mid_cloud',
        'lcdclcll':'low_cloud',
        'vissfc':'sfcvis',
        'hgt263_k':'hgt_m10c',
        'hgt253_k':'hgt_m20c',
        'ltngclm':'lightning',
        'sbt124toa':'simsat',
        'hgt0c':'0chgt'
    })
    zH5 = data['temperature'].squeeze()
    zH5_crs = zH5.metpy.cartopy_crs

    vis = data['sfcvis'].squeeze()*0.000621371
    reflectivity = data['radar'].squeeze()
    cape = data['cape'].squeeze()
    lightning=data['lightning'].squeeze()
    dgz_depth = data['hgt_m20c'].squeeze()-data['hgt_m10c'].squeeze()
    simsat = data['simsat'].squeeze()
    hgt0c = data['0chgt'].squeeze()*3.28084
    print(np.max(simsat))
    print(np.min(simsat))
    print(lightning)
    print(np.max(lightning))

    vertical, = data['temperature'].metpy.coordinates('vertical')
    time = data['temperature'].metpy.time
    x, y = data['temperature'].metpy.coordinates('x', 'y')
    lat, lon = xr.broadcast(y, x)

    fig = plt.figure(figsize=(42,15))

    gs = fig.add_gridspec(ncols=3,nrows= 2, width_ratios=[1,2,1])
    ax1 = fig.add_subplot(gs[:, 1], projection = zH5_crs)
    ax2 = fig.add_subplot(gs[0, 0], projection = zH5_crs)
    ax3 = fig.add_subplot(gs[1, 0], projection = zH5_crs)
    ax4 = fig.add_subplot(gs[0, 2], projection = zH5_crs)
    ax5 = fig.add_subplot(gs[1, 2], projection = zH5_crs)

    ax1.coastlines(resolution='10m')
    ax1.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax1.add_feature(cfeature.STATES.with_scale('10m'))

    ax2.coastlines(resolution='10m')
    ax2.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax2.add_feature(cfeature.STATES.with_scale('10m'))

    ax3.coastlines(resolution='10m')
    ax3.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax3.add_feature(cfeature.STATES.with_scale('10m'))

    ax4.coastlines(resolution='10m')
    ax4.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax4.add_feature(cfeature.STATES.with_scale('10m'))

    ax5.coastlines(resolution='10m')
    ax5.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax5.add_feature(cfeature.STATES.with_scale('10m'))

    #fig.suptitle("HRRR Forecast valid at " + time[0].dt.strftime('%Y-%m-%d %H:%MZ').item(),fontsize=36)

    dtfs = str(time.dt.strftime('%Y-%m-%d_%H%MZ').item())
    print(dtfs)

    vis = ax1.contourf(x,y,vis,levels=[0.00625,0.125,0.25,0.5,0.75,1,2,3,4,5],extend='min',colors=['#ffb3ff','magenta','deeppink','hotpink','mediumvioletred','crimson','orangered','darkorange','orange','gold'],alpha=0.7)
    cbr = fig.colorbar(vis,orientation = 'horizontal', aspect = 80, ax = ax1, pad = 0.01,
                        extendrect=False, ticks = [0.00625,0.125,0.25,0.5,0.75,1,2,3,4,5], shrink=0.7)
    cbr.set_label('Surface Visibility (mi)')

    refp = ax4.contourf(x,y,reflectivity, levels=[20, 25, 30, 35, 40, 45, 50, 55, 60, 65], alpha = 0.7, cmap = 'Greens',transform=zH5_crs) #colors=['#0099ff00', '#4D8080ff', '#666666ff', '#804d4dff','#993333ff','#B33333ff','#CC1a1aff','#E60000ff','#0000e6','#0000cc','#0000b3','#2d00b3','#5900b3','#8600b3','#b300b3','#b30086'])
    capep = ax4.contourf(x, y, cape, levels=[100, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2500, 3000], alpha = 0.7, cmap='RdPu')#['#0099ff00', '#4066ffb3', '#8066ff8c', '#BF66ff66','#8cff66','#b3ff66','#d9ff66','#ffff66','#ffd966','#ffcc66','#ffb366','#ff8c66','#ff6666','#ff668c','#ff66b3','#ff66d9','#ff66ff'])
    lgt = ax4.contour(x,y,lightning,levels=[0.5,1,1.5,2,2.5,3,3.5,4,4.5,5])
    cb = fig.colorbar(capep, orientation='vertical', pad = 0.01, aspect = 50, ax = ax4, extendrect=False, ticks=[100, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2500, 3000])
    cb.set_label('CAPE (J/kg)', size='large')

    dgz = ax5.contourf(x,y,dgz_depth,alpha=0.7)
    cbar = fig.colorbar(dgz,orientation='vertical', pad = 0.01, aspect = 50, ax = ax5, extendrect=False)

    ir = ax2.contourf(x,y,simsat,cmap='gnuplot2',levs=range(173,263,5),extend='both')
    cbar = fig.colorbar(ir,orientation='vertical', pad = 0.01, aspect = 50, ax = ax2, extendrect=False)

    hgt_frz = ax3.contourf(x,y,hgt0c,cmap='rainbow',levels=range(0,10000,500),extend='max')
    cbar2 = fig.colorbar(hgt_frz,orientation='vertical', pad = 0.01, aspect = 50, ax = ax3, extendrect=False)

    sub_w1 = 260
    sub_w = 262
    sub_e = 295
    sub_n = 50
    sub_s = 25

    ax1.set_extent((sub_w1, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax2.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax3.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax4.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax5.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot

    fig.tight_layout()
    plt.savefig(output_dir+'/HRRR_test/newparams6_'+dtfs+'.png')

847f9c973fde475d3fb2847452e8f206286ce48f
<mask token>


def mkdir_p(mypath):
    """Creates a directory. equivalent to using mkdir -p on the command line"""
    from errno import EEXIST
    from os import makedirs, path
    try:
        makedirs(mypath)
    except OSError as exc:
        if exc.errno == EEXIST and path.isdir(mypath):
            pass
        else:
            raise


<mask token>


def get_init_hr(hour):
    if int(hour) < 3:
        init_hour = '00'
    elif int(hour) < 9:
        init_hour = '06'
    elif int(hour) < 15:
        init_hour = '12'
    elif int(hour) < 21:
        init_hour = '18'
    else:
        init_hour = '00'
    return init_hour


<mask token>

<mask token>


def mkdir_p(mypath):
    """Creates a directory. equivalent to using mkdir -p on the command line"""
    from errno import EEXIST
    from os import makedirs, path
    try:
        makedirs(mypath)
    except OSError as exc:
        if exc.errno == EEXIST and path.isdir(mypath):
            pass
        else:
            raise


<mask token>
if startTime.month < 10:
    month = '0' + str(startTime.month)
else:
    month = str(startTime.month)
if startTime.day < 10:
    day = '0' + str(startTime.day)
else:
    day = str(startTime.day)
if startTime.hour < 10:
    hour = '0' + str(startTime.hour)
else:
    hour = str(startTime.hour)
<mask token>


def get_init_hr(hour):
    if int(hour) < 3:
        init_hour = '00'
    elif int(hour) < 9:
        init_hour = '06'
    elif int(hour) < 15:
        init_hour = '12'
    elif int(hour) < 21:
        init_hour = '18'
    else:
        init_hour = '00'
    return init_hour


<mask token>
print(url)
<mask token>
mkdir_p(output_dir)
mkdir_p(output_dir + '/HRRR_test')
<mask token>
for i in range(1, 49):
    fc_hr = init_hr + dt.timedelta(hours=1 * i)
    forecast_hour = times[0].values
    data = ds.metpy.parse_cf()
    data = data.isel(time=i)
    data = data.rename({'cfrzrsfc': 'catice', 'cicepsfc': 'catsleet',
        'crainsfc': 'catrain', 'csnowsfc': 'catsnow', 'tcdcclm': 'tcc',
        'tmpprs': 'temperature', 'ugrd10m': 'u', 'vgrd10m': 'v', 'mslmamsl':
        'mslp', 'tmp2m': 'sfc_temp', 'dpt2m': 'sfc_td', 'refcclm': 'radar',
        'apcpsfc': 'qpf', 'capesfc': 'cape', 'gustsfc': 'sfcgust',
        'hcdchcll': 'high_cloud', 'mcdcmcll': 'mid_cloud', 'lcdclcll':
        'low_cloud', 'vissfc': 'sfcvis', 'hgt263_k': 'hgt_m10c', 'hgt253_k':
        'hgt_m20c', 'ltngclm': 'lightning', 'sbt124toa': 'simsat', 'hgt0c':
        '0chgt'})
    zH5 = data['temperature'].squeeze()
    zH5_crs = zH5.metpy.cartopy_crs
    vis = data['sfcvis'].squeeze() * 0.000621371
    reflectivity = data['radar'].squeeze()
    cape = data['cape'].squeeze()
    lightning = data['lightning'].squeeze()
    dgz_depth = data['hgt_m20c'].squeeze() - data['hgt_m10c'].squeeze()
    simsat = data['simsat'].squeeze()
    hgt0c = data['0chgt'].squeeze() * 3.28084
    print(np.max(simsat))
    print(np.min(simsat))
    print(lightning)
    print(np.max(lightning))
    vertical, = data['temperature'].metpy.coordinates('vertical')
    time = data['temperature'].metpy.time
    x, y = data['temperature'].metpy.coordinates('x', 'y')
    lat, lon = xr.broadcast(y, x)
    fig = plt.figure(figsize=(42, 15))
    gs = fig.add_gridspec(ncols=3, nrows=2, width_ratios=[1, 2, 1])
    ax1 = fig.add_subplot(gs[:, 1], projection=zH5_crs)
    ax2 = fig.add_subplot(gs[0, 0], projection=zH5_crs)
    ax3 = fig.add_subplot(gs[1, 0], projection=zH5_crs)
    ax4 = fig.add_subplot(gs[0, 2], projection=zH5_crs)
    ax5 = fig.add_subplot(gs[1, 2], projection=zH5_crs)
    ax1.coastlines(resolution='10m')
    ax1.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax1.add_feature(cfeature.STATES.with_scale('10m'))
    ax2.coastlines(resolution='10m')
    ax2.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax2.add_feature(cfeature.STATES.with_scale('10m'))
    ax3.coastlines(resolution='10m')
    ax3.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax3.add_feature(cfeature.STATES.with_scale('10m'))
    ax4.coastlines(resolution='10m')
    ax4.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax4.add_feature(cfeature.STATES.with_scale('10m'))
    ax5.coastlines(resolution='10m')
    ax5.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax5.add_feature(cfeature.STATES.with_scale('10m'))
    dtfs = str(time.dt.strftime('%Y-%m-%d_%H%MZ').item())
    print(dtfs)
    vis = ax1.contourf(x, y, vis, levels=[0.00625, 0.125, 0.25, 0.5, 0.75, 
        1, 2, 3, 4, 5], extend='min', colors=['#ffb3ff', 'magenta',
        'deeppink', 'hotpink', 'mediumvioletred', 'crimson', 'orangered',
        'darkorange', 'orange', 'gold'], alpha=0.7)
    cbr = fig.colorbar(vis, orientation='horizontal', aspect=80, ax=ax1,
        pad=0.01, extendrect=False, ticks=[0.00625, 0.125, 0.25, 0.5, 0.75,
        1, 2, 3, 4, 5], shrink=0.7)
    cbr.set_label('Surface Visibility (mi)')
    refp = ax4.contourf(x, y, reflectivity, levels=[20, 25, 30, 35, 40, 45,
        50, 55, 60, 65], alpha=0.7, cmap='Greens', transform=zH5_crs)
    capep = ax4.contourf(x, y, cape, levels=[100, 250, 500, 750, 1000, 1250,
        1500, 1750, 2000, 2500, 3000], alpha=0.7, cmap='RdPu')
    lgt = ax4.contour(x, y, lightning, levels=[0.5, 1, 1.5, 2, 2.5, 3, 3.5,
        4, 4.5, 5])
    cb = fig.colorbar(capep, orientation='vertical', pad=0.01, aspect=50,
        ax=ax4, extendrect=False, ticks=[100, 250, 500, 750, 1000, 1250, 
        1500, 1750, 2000, 2500, 3000])
    cb.set_label('CAPE (J/kg)', size='large')
    dgz = ax5.contourf(x, y, dgz_depth, alpha=0.7)
    cbar = fig.colorbar(dgz, orientation='vertical', pad=0.01, aspect=50,
        ax=ax5, extendrect=False)
    ir = ax2.contourf(x, y, simsat, cmap='gnuplot2', levs=range(173, 263, 5
        ), extend='both')
    cbar = fig.colorbar(ir, orientation='vertical', pad=0.01, aspect=50, ax
        =ax2, extendrect=False)
    hgt_frz = ax3.contourf(x, y, hgt0c, cmap='rainbow', levels=range(0, 
        10000, 500), extend='max')
    cbar2 = fig.colorbar(hgt_frz, orientation='vertical', pad=0.01, aspect=
        50, ax=ax3, extendrect=False)
    sub_w1 = 260
    sub_w = 262
    sub_e = 295
    sub_n = 50
    sub_s = 25
    ax1.set_extent((sub_w1, sub_e, sub_s, sub_n))
    ax2.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax3.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax4.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax5.set_extent((sub_w, sub_e, sub_s, sub_n))
    fig.tight_layout()
    plt.savefig(output_dir + '/HRRR_test/newparams6_' + dtfs + '.png')

<mask token>


def mkdir_p(mypath):
    """Creates a directory. equivalent to using mkdir -p on the command line"""
    from errno import EEXIST
    from os import makedirs, path
    try:
        makedirs(mypath)
    except OSError as exc:
        if exc.errno == EEXIST and path.isdir(mypath):
            pass
        else:
            raise


startTime = datetime.now()
year = startTime.year
if startTime.month < 10:
    month = '0' + str(startTime.month)
else:
    month = str(startTime.month)
if startTime.day < 10:
    day = '0' + str(startTime.day)
else:
    day = str(startTime.day)
if startTime.hour < 10:
    hour = '0' + str(startTime.hour)
else:
    hour = str(startTime.hour)
mdate = str(year) + str(month) + str(day)


def get_init_hr(hour):
    if int(hour) < 3:
        init_hour = '00'
    elif int(hour) < 9:
        init_hour = '06'
    elif int(hour) < 15:
        init_hour = '12'
    elif int(hour) < 21:
        init_hour = '18'
    else:
        init_hour = '00'
    return init_hour


init_hour = get_init_hr(hour)
url = ('http://nomads.ncep.noaa.gov:80/dods/hrrr/hrrr' + mdate +
    '/hrrr_sfc.t' + get_init_hr(hour) + 'z')
print(url)
output_dir = str(year) + str(month) + str(day) + '_' + str(init_hour) + '00'
mkdir_p(output_dir)
mkdir_p(output_dir + '/HRRR_test')
ds = xr.open_dataset(url)
init_hr = dt.datetime(int(year), int(month), int(day), int(init_hour))
times = ds['tmp2m'].metpy.time
init_time = ds['time'][0]
lats = np.arange(25, 55, 0.25)
lons = np.arange(260, 310, 0.25)
for i in range(1, 49):
    fc_hr = init_hr + dt.timedelta(hours=1 * i)
    forecast_hour = times[0].values
    data = ds.metpy.parse_cf()
    data = data.isel(time=i)
    data = data.rename({'cfrzrsfc': 'catice', 'cicepsfc': 'catsleet',
        'crainsfc': 'catrain', 'csnowsfc': 'catsnow', 'tcdcclm': 'tcc',
        'tmpprs': 'temperature', 'ugrd10m': 'u', 'vgrd10m': 'v', 'mslmamsl':
        'mslp', 'tmp2m': 'sfc_temp', 'dpt2m': 'sfc_td', 'refcclm': 'radar',
        'apcpsfc': 'qpf', 'capesfc': 'cape', 'gustsfc': 'sfcgust',
        'hcdchcll': 'high_cloud', 'mcdcmcll': 'mid_cloud', 'lcdclcll':
        'low_cloud', 'vissfc': 'sfcvis', 'hgt263_k': 'hgt_m10c', 'hgt253_k':
        'hgt_m20c', 'ltngclm': 'lightning', 'sbt124toa': 'simsat', 'hgt0c':
        '0chgt'})
    zH5 = data['temperature'].squeeze()
    zH5_crs = zH5.metpy.cartopy_crs
    vis = data['sfcvis'].squeeze() * 0.000621371
    reflectivity = data['radar'].squeeze()
    cape = data['cape'].squeeze()
    lightning = data['lightning'].squeeze()
    dgz_depth = data['hgt_m20c'].squeeze() - data['hgt_m10c'].squeeze()
    simsat = data['simsat'].squeeze()
    hgt0c = data['0chgt'].squeeze() * 3.28084
    print(np.max(simsat))
    print(np.min(simsat))
    print(lightning)
    print(np.max(lightning))
    vertical, = data['temperature'].metpy.coordinates('vertical')
    time = data['temperature'].metpy.time
    x, y = data['temperature'].metpy.coordinates('x', 'y')
    lat, lon = xr.broadcast(y, x)
    fig = plt.figure(figsize=(42, 15))
    gs = fig.add_gridspec(ncols=3, nrows=2, width_ratios=[1, 2, 1])
    ax1 = fig.add_subplot(gs[:, 1], projection=zH5_crs)
    ax2 = fig.add_subplot(gs[0, 0], projection=zH5_crs)
    ax3 = fig.add_subplot(gs[1, 0], projection=zH5_crs)
    ax4 = fig.add_subplot(gs[0, 2], projection=zH5_crs)
    ax5 = fig.add_subplot(gs[1, 2], projection=zH5_crs)
    ax1.coastlines(resolution='10m')
    ax1.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax1.add_feature(cfeature.STATES.with_scale('10m'))
    ax2.coastlines(resolution='10m')
    ax2.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax2.add_feature(cfeature.STATES.with_scale('10m'))
    ax3.coastlines(resolution='10m')
    ax3.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax3.add_feature(cfeature.STATES.with_scale('10m'))
    ax4.coastlines(resolution='10m')
    ax4.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax4.add_feature(cfeature.STATES.with_scale('10m'))
    ax5.coastlines(resolution='10m')
    ax5.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax5.add_feature(cfeature.STATES.with_scale('10m'))
    dtfs = str(time.dt.strftime('%Y-%m-%d_%H%MZ').item())
    print(dtfs)
    vis = ax1.contourf(x, y, vis, levels=[0.00625, 0.125, 0.25, 0.5, 0.75, 
        1, 2, 3, 4, 5], extend='min', colors=['#ffb3ff', 'magenta',
        'deeppink', 'hotpink', 'mediumvioletred', 'crimson', 'orangered',
        'darkorange', 'orange', 'gold'], alpha=0.7)
    cbr = fig.colorbar(vis, orientation='horizontal', aspect=80, ax=ax1,
        pad=0.01, extendrect=False, ticks=[0.00625, 0.125, 0.25, 0.5, 0.75,
        1, 2, 3, 4, 5], shrink=0.7)
    cbr.set_label('Surface Visibility (mi)')
    refp = ax4.contourf(x, y, reflectivity, levels=[20, 25, 30, 35, 40, 45,
        50, 55, 60, 65], alpha=0.7, cmap='Greens', transform=zH5_crs)
    capep = ax4.contourf(x, y, cape, levels=[100, 250, 500, 750, 1000, 1250,
        1500, 1750, 2000, 2500, 3000], alpha=0.7, cmap='RdPu')
    lgt = ax4.contour(x, y, lightning, levels=[0.5, 1, 1.5, 2, 2.5, 3, 3.5,
        4, 4.5, 5])
    cb = fig.colorbar(capep, orientation='vertical', pad=0.01, aspect=50,
        ax=ax4, extendrect=False, ticks=[100, 250, 500, 750, 1000, 1250, 
        1500, 1750, 2000, 2500, 3000])
    cb.set_label('CAPE (J/kg)', size='large')
    dgz = ax5.contourf(x, y, dgz_depth, alpha=0.7)
    cbar = fig.colorbar(dgz, orientation='vertical', pad=0.01, aspect=50,
        ax=ax5, extendrect=False)
    ir = ax2.contourf(x, y, simsat, cmap='gnuplot2', levs=range(173, 263, 5
        ), extend='both')
    cbar = fig.colorbar(ir, orientation='vertical', pad=0.01, aspect=50, ax
        =ax2, extendrect=False)
    hgt_frz = ax3.contourf(x, y, hgt0c, cmap='rainbow', levels=range(0, 
        10000, 500), extend='max')
    cbar2 = fig.colorbar(hgt_frz, orientation='vertical', pad=0.01, aspect=
        50, ax=ax3, extendrect=False)
    sub_w1 = 260
    sub_w = 262
    sub_e = 295
    sub_n = 50
    sub_s = 25
    ax1.set_extent((sub_w1, sub_e, sub_s, sub_n))
    ax2.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax3.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax4.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax5.set_extent((sub_w, sub_e, sub_s, sub_n))
    fig.tight_layout()
    plt.savefig(output_dir + '/HRRR_test/newparams6_' + dtfs + '.png')

import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from netCDF4 import num2date
import numpy as np
import xarray as xr
from siphon.catalog import TDSCatalog
from datetime import datetime
import datetime as dt
from xarray.backends import NetCDF4DataStore
import cartopy
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.ndimage import gaussian_filter
import metpy.calc as mpcalc
import numpy.ma as ma
from metpy.units import units
import scipy.ndimage as ndimage
from metpy.plots import USCOUNTIES
import matplotlib.patches as mpatches


def mkdir_p(mypath):
    """Creates a directory. equivalent to using mkdir -p on the command line"""
    from errno import EEXIST
    from os import makedirs, path
    try:
        makedirs(mypath)
    except OSError as exc:
        if exc.errno == EEXIST and path.isdir(mypath):
            pass
        else:
            raise


startTime = datetime.now()
year = startTime.year
if startTime.month < 10:
    month = '0' + str(startTime.month)
else:
    month = str(startTime.month)
if startTime.day < 10:
    day = '0' + str(startTime.day)
else:
    day = str(startTime.day)
if startTime.hour < 10:
    hour = '0' + str(startTime.hour)
else:
    hour = str(startTime.hour)
mdate = str(year) + str(month) + str(day)


def get_init_hr(hour):
    if int(hour) < 3:
        init_hour = '00'
    elif int(hour) < 9:
        init_hour = '06'
    elif int(hour) < 15:
        init_hour = '12'
    elif int(hour) < 21:
        init_hour = '18'
    else:
        init_hour = '00'
    return init_hour


init_hour = get_init_hr(hour)
url = ('http://nomads.ncep.noaa.gov:80/dods/hrrr/hrrr' + mdate +
    '/hrrr_sfc.t' + get_init_hr(hour) + 'z')
print(url)
output_dir = str(year) + str(month) + str(day) + '_' + str(init_hour) + '00'
mkdir_p(output_dir)
mkdir_p(output_dir + '/HRRR_test')
ds = xr.open_dataset(url)
init_hr = dt.datetime(int(year), int(month), int(day), int(init_hour))
times = ds['tmp2m'].metpy.time
init_time = ds['time'][0]
lats = np.arange(25, 55, 0.25)
lons = np.arange(260, 310, 0.25)
for i in range(1, 49):
    fc_hr = init_hr + dt.timedelta(hours=1 * i)
    forecast_hour = times[0].values
    data = ds.metpy.parse_cf()
    data = data.isel(time=i)
    data = data.rename({'cfrzrsfc': 'catice', 'cicepsfc': 'catsleet',
        'crainsfc': 'catrain', 'csnowsfc': 'catsnow', 'tcdcclm': 'tcc',
        'tmpprs': 'temperature', 'ugrd10m': 'u', 'vgrd10m': 'v', 'mslmamsl':
        'mslp', 'tmp2m': 'sfc_temp', 'dpt2m': 'sfc_td', 'refcclm': 'radar',
        'apcpsfc': 'qpf', 'capesfc': 'cape', 'gustsfc': 'sfcgust',
        'hcdchcll': 'high_cloud', 'mcdcmcll': 'mid_cloud', 'lcdclcll':
        'low_cloud', 'vissfc': 'sfcvis', 'hgt263_k': 'hgt_m10c', 'hgt253_k':
        'hgt_m20c', 'ltngclm': 'lightning', 'sbt124toa': 'simsat', 'hgt0c':
        '0chgt'})
    zH5 = data['temperature'].squeeze()
    zH5_crs = zH5.metpy.cartopy_crs
    vis = data['sfcvis'].squeeze() * 0.000621371
    reflectivity = data['radar'].squeeze()
    cape = data['cape'].squeeze()
    lightning = data['lightning'].squeeze()
    dgz_depth = data['hgt_m20c'].squeeze() - data['hgt_m10c'].squeeze()
    simsat = data['simsat'].squeeze()
    hgt0c = data['0chgt'].squeeze() * 3.28084
    print(np.max(simsat))
    print(np.min(simsat))
    print(lightning)
    print(np.max(lightning))
    vertical, = data['temperature'].metpy.coordinates('vertical')
    time = data['temperature'].metpy.time
    x, y = data['temperature'].metpy.coordinates('x', 'y')
    lat, lon = xr.broadcast(y, x)
    fig = plt.figure(figsize=(42, 15))
    gs = fig.add_gridspec(ncols=3, nrows=2, width_ratios=[1, 2, 1])
    ax1 = fig.add_subplot(gs[:, 1], projection=zH5_crs)
    ax2 = fig.add_subplot(gs[0, 0], projection=zH5_crs)
    ax3 = fig.add_subplot(gs[1, 0], projection=zH5_crs)
    ax4 = fig.add_subplot(gs[0, 2], projection=zH5_crs)
    ax5 = fig.add_subplot(gs[1, 2], projection=zH5_crs)
    ax1.coastlines(resolution='10m')
    ax1.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax1.add_feature(cfeature.STATES.with_scale('10m'))
    ax2.coastlines(resolution='10m')
    ax2.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax2.add_feature(cfeature.STATES.with_scale('10m'))
    ax3.coastlines(resolution='10m')
    ax3.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax3.add_feature(cfeature.STATES.with_scale('10m'))
    ax4.coastlines(resolution='10m')
    ax4.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax4.add_feature(cfeature.STATES.with_scale('10m'))
    ax5.coastlines(resolution='10m')
    ax5.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax5.add_feature(cfeature.STATES.with_scale('10m'))
    dtfs = str(time.dt.strftime('%Y-%m-%d_%H%MZ').item())
    print(dtfs)
    vis = ax1.contourf(x, y, vis, levels=[0.00625, 0.125, 0.25, 0.5, 0.75, 
        1, 2, 3, 4, 5], extend='min', colors=['#ffb3ff', 'magenta',
        'deeppink', 'hotpink', 'mediumvioletred', 'crimson', 'orangered',
        'darkorange', 'orange', 'gold'], alpha=0.7)
    cbr = fig.colorbar(vis, orientation='horizontal', aspect=80, ax=ax1,
        pad=0.01, extendrect=False, ticks=[0.00625, 0.125, 0.25, 0.5, 0.75,
        1, 2, 3, 4, 5], shrink=0.7)
    cbr.set_label('Surface Visibility (mi)')
    refp = ax4.contourf(x, y, reflectivity, levels=[20, 25, 30, 35, 40, 45,
        50, 55, 60, 65], alpha=0.7, cmap='Greens', transform=zH5_crs)
    capep = ax4.contourf(x, y, cape, levels=[100, 250, 500, 750, 1000, 1250,
        1500, 1750, 2000, 2500, 3000], alpha=0.7, cmap='RdPu')
    lgt = ax4.contour(x, y, lightning, levels=[0.5, 1, 1.5, 2, 2.5, 3, 3.5,
        4, 4.5, 5])
    cb = fig.colorbar(capep, orientation='vertical', pad=0.01, aspect=50,
        ax=ax4, extendrect=False, ticks=[100, 250, 500, 750, 1000, 1250, 
        1500, 1750, 2000, 2500, 3000])
    cb.set_label('CAPE (J/kg)', size='large')
    dgz = ax5.contourf(x, y, dgz_depth, alpha=0.7)
    cbar = fig.colorbar(dgz, orientation='vertical', pad=0.01, aspect=50,
        ax=ax5, extendrect=False)
    ir = ax2.contourf(x, y, simsat, cmap='gnuplot2', levs=range(173, 263, 5
        ), extend='both')
    cbar = fig.colorbar(ir, orientation='vertical', pad=0.01, aspect=50, ax
        =ax2, extendrect=False)
    hgt_frz = ax3.contourf(x, y, hgt0c, cmap='rainbow', levels=range(0, 
        10000, 500), extend='max')
    cbar2 = fig.colorbar(hgt_frz, orientation='vertical', pad=0.01, aspect=
        50, ax=ax3, extendrect=False)
    sub_w1 = 260
    sub_w = 262
    sub_e = 295
    sub_n = 50
    sub_s = 25
    ax1.set_extent((sub_w1, sub_e, sub_s, sub_n))
    ax2.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax3.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax4.set_extent((sub_w, sub_e, sub_s, sub_n))
    ax5.set_extent((sub_w, sub_e, sub_s, sub_n))
    fig.tight_layout()
    plt.savefig(output_dir + '/HRRR_test/newparams6_' + dtfs + '.png')

import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from netCDF4 import num2date
import numpy as np
import xarray as xr
from siphon.catalog import TDSCatalog
from datetime import datetime
import datetime as dt
from xarray.backends import NetCDF4DataStore
import cartopy
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.ndimage import gaussian_filter
import metpy.calc as mpcalc
import numpy.ma as ma
from metpy.units import units
import scipy.ndimage as ndimage
from metpy.plots import USCOUNTIES
import matplotlib.patches as mpatches

def mkdir_p(mypath):
    '''Creates a directory. equivalent to using mkdir -p on the command line'''

    from errno import EEXIST
    from os import makedirs,path

    try:
        makedirs(mypath)
    except OSError as exc: # Python >2.5
        if exc.errno == EEXIST and path.isdir(mypath):
            pass
        else: raise

startTime=datetime.now()

year = startTime.year

if startTime.month <10:
    month = '0'+str(startTime.month)
else:
    month = str(startTime.month)

if startTime.day <10:
    day = '0'+str(startTime.day)
else:
    day = str(startTime.day)

if startTime.hour <10:
    hour = '0'+str(startTime.hour)
else:
    hour = str(startTime.hour)

mdate = str(year)+str(month)+str(day)

def get_init_hr(hour):
    if int(hour) <3:
        init_hour = '00'
    elif int(hour) <9:
        init_hour = '06'
    elif int(hour) <15:
        init_hour = '12'
    elif int(hour) <21:
        init_hour = '18'
    else:
        init_hour = '00'
    return(init_hour)
init_hour = get_init_hr(hour)
url = 'http://nomads.ncep.noaa.gov:80/dods/hrrr/hrrr'+mdate+'/hrrr_sfc.t'+get_init_hr(hour)+'z'
#url='http://nomads.ncep.noaa.gov:80/dods/hrrr/hrrr20201231/hrrr_sfc.t00z'
print(url)

# Create new directory
output_dir = str(year)+str(month)+str(day)+'_'+str(init_hour)+'00'
#output_dir = '20201231_0000'
mkdir_p(output_dir)
mkdir_p(output_dir+'/HRRR_test')
#Parse data using MetPy
ds = xr.open_dataset(url)
init_hr = dt.datetime(int(year),int(month),int(day),int(init_hour))
times = ds['tmp2m'].metpy.time
init_time = ds['time'][0]

lats = np.arange(25,55,0.25)
lons = np.arange(260,310,0.25)

for i in range(1,49):
    fc_hr = init_hr+dt.timedelta(hours=1*i)
    forecast_hour = times[0].values

    data = ds.metpy.parse_cf()
    data = data.isel(time=i)

    #Rename variables to useful things
    data = data.rename({
        'cfrzrsfc':'catice',
        'cicepsfc':'catsleet',
        'crainsfc':'catrain',
        'csnowsfc':'catsnow',
        'tcdcclm':'tcc',
        'tmpprs': 'temperature',
        'ugrd10m': 'u',
        'vgrd10m': 'v',
        'mslmamsl':'mslp',
        'tmp2m':'sfc_temp',
        'dpt2m':'sfc_td',
        'refcclm':'radar',
        'apcpsfc':'qpf',
        'capesfc':'cape',
        'gustsfc':'sfcgust',
        'hcdchcll':'high_cloud',
        'mcdcmcll':'mid_cloud',
        'lcdclcll':'low_cloud',
        'vissfc':'sfcvis',
        'hgt263_k':'hgt_m10c',
        'hgt253_k':'hgt_m20c',
        'ltngclm':'lightning',
        'sbt124toa':'simsat',
        'hgt0c':'0chgt'
    })
    zH5 = data['temperature'].squeeze()
    zH5_crs = zH5.metpy.cartopy_crs

    vis = data['sfcvis'].squeeze()*0.000621371
    reflectivity = data['radar'].squeeze()
    cape = data['cape'].squeeze()
    lightning=data['lightning'].squeeze()
    dgz_depth = data['hgt_m20c'].squeeze()-data['hgt_m10c'].squeeze()
    simsat = data['simsat'].squeeze()
    hgt0c = data['0chgt'].squeeze()*3.28084
    print(np.max(simsat))
    print(np.min(simsat))
    print(lightning)
    print(np.max(lightning))

    vertical, = data['temperature'].metpy.coordinates('vertical')
    time = data['temperature'].metpy.time
    x, y = data['temperature'].metpy.coordinates('x', 'y')
    lat, lon = xr.broadcast(y, x)

    fig = plt.figure(figsize=(42,15))

    gs = fig.add_gridspec(ncols=3,nrows= 2, width_ratios=[1,2,1])
    ax1 = fig.add_subplot(gs[:, 1], projection = zH5_crs)
    ax2 = fig.add_subplot(gs[0, 0], projection = zH5_crs)
    ax3 = fig.add_subplot(gs[1, 0], projection = zH5_crs)
    ax4 = fig.add_subplot(gs[0, 2], projection = zH5_crs)
    ax5 = fig.add_subplot(gs[1, 2], projection = zH5_crs)

    ax1.coastlines(resolution='10m')
    ax1.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax1.add_feature(cfeature.STATES.with_scale('10m'))

    ax2.coastlines(resolution='10m')
    ax2.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax2.add_feature(cfeature.STATES.with_scale('10m'))

    ax3.coastlines(resolution='10m')
    ax3.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax3.add_feature(cfeature.STATES.with_scale('10m'))

    ax4.coastlines(resolution='10m')
    ax4.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax4.add_feature(cfeature.STATES.with_scale('10m'))

    ax5.coastlines(resolution='10m')
    ax5.add_feature(cfeature.BORDERS.with_scale('10m'))
    ax5.add_feature(cfeature.STATES.with_scale('10m'))

    #fig.suptitle("HRRR Forecast valid at " + time[0].dt.strftime('%Y-%m-%d %H:%MZ').item(),fontsize=36)

    dtfs = str(time.dt.strftime('%Y-%m-%d_%H%MZ').item())
    print(dtfs)

    vis = ax1.contourf(x,y,vis,levels=[0.00625,0.125,0.25,0.5,0.75,1,2,3,4,5],extend='min',colors=['#ffb3ff','magenta','deeppink','hotpink','mediumvioletred','crimson','orangered','darkorange','orange','gold'],alpha=0.7)
    cbr = fig.colorbar(vis,orientation = 'horizontal', aspect = 80, ax = ax1, pad = 0.01,
                        extendrect=False, ticks = [0.00625,0.125,0.25,0.5,0.75,1,2,3,4,5], shrink=0.7)
    cbr.set_label('Surface Visibility (mi)')

    refp = ax4.contourf(x,y,reflectivity, levels=[20, 25, 30, 35, 40, 45, 50, 55, 60, 65], alpha = 0.7, cmap = 'Greens',transform=zH5_crs) #colors=['#0099ff00', '#4D8080ff', '#666666ff', '#804d4dff','#993333ff','#B33333ff','#CC1a1aff','#E60000ff','#0000e6','#0000cc','#0000b3','#2d00b3','#5900b3','#8600b3','#b300b3','#b30086'])
    capep = ax4.contourf(x, y, cape, levels=[100, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2500, 3000], alpha = 0.7, cmap='RdPu')#['#0099ff00', '#4066ffb3', '#8066ff8c', '#BF66ff66','#8cff66','#b3ff66','#d9ff66','#ffff66','#ffd966','#ffcc66','#ffb366','#ff8c66','#ff6666','#ff668c','#ff66b3','#ff66d9','#ff66ff'])
    lgt = ax4.contour(x,y,lightning,levels=[0.5,1,1.5,2,2.5,3,3.5,4,4.5,5])
    cb = fig.colorbar(capep, orientation='vertical', pad = 0.01, aspect = 50, ax = ax4, extendrect=False, ticks=[100, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2500, 3000])
    cb.set_label('CAPE (J/kg)', size='large')

    dgz = ax5.contourf(x,y,dgz_depth,alpha=0.7)
    cbar = fig.colorbar(dgz,orientation='vertical', pad = 0.01, aspect = 50, ax = ax5, extendrect=False)

    ir = ax2.contourf(x,y,simsat,cmap='gnuplot2',levs=range(173,263,5),extend='both')
    cbar = fig.colorbar(ir,orientation='vertical', pad = 0.01, aspect = 50, ax = ax2, extendrect=False)

    hgt_frz = ax3.contourf(x,y,hgt0c,cmap='rainbow',levels=range(0,10000,500),extend='max')
    cbar2 = fig.colorbar(hgt_frz,orientation='vertical', pad = 0.01, aspect = 50, ax = ax3, extendrect=False)

    sub_w1 = 260
    sub_w = 262
    sub_e = 295
    sub_n = 50
    sub_s = 25

    ax1.set_extent((sub_w1, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax2.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax3.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax4.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot
    ax5.set_extent((sub_w, sub_e, sub_s, sub_n))#, crs = zH5_crs)    # Set a title and show the plot

    fig.tight_layout()
    plt.savefig(output_dir+'/HRRR_test/newparams6_'+dtfs+'.png')
