initial commit

input_shapefile/NCCtyMerc.cpg

+UTF-8
\ No newline at end of file

input_shapefile/NCCtyMerc.prj

+PROJCS["WGS_1984_World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]
\ No newline at end of file

input_shapefile/NCSteMerc.cpg

+UTF-8
\ No newline at end of file

input_shapefile/NCSteMerc.prj

+PROJCS["WGS_1984_World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]
\ No newline at end of file

input_shapefile/NCSteMerc.shp.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<metadata xml:lang="en"><Esri><CreaDate>20210818</CreaDate><CreaTime>10241200</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>TRUE</SyncOnce><DataProperties><lineage><Process ToolSource="c:\program files (x86)\arcgis\desktop10.8\ArcToolbox\Toolboxes\Data Management Tools.tbx\Project" Date="20210818" Time="102412">Project Export_Output E:\jared\Shapefiles\NCSteMerc.shp PROJCS['WGS_1984_World_Mercator',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],UNIT['Meter',1.0]] WGS_1984_(ITRF00)_To_NAD_1983 GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]] NO_PRESERVE_SHAPE # NO_VERTICAL</Process></lineage></DataProperties></Esri></metadata>

plot_hrrr_nc-gridded.py

+################################################################################
+# Written by Jared Rennie, NCICS
+# Derived this iPython Notebook, from Kevin Tyle
+# https://github.com/ktyle/python_pangeo_ams2021/blob/main/notebooks/03_Pangeo_HRRR.ipynb
+################################################################################
+
+# Import Packages
+import sys, time, datetime, os
+import numpy as np
+import numpy.ma as ma
+import geopandas as gpd
+import pandas as pd
+import xarray as xr
+import s3fs
+import metpy
+
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+
+import cartopy
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import cartopy.io.shapereader as shpreader
+
+import rioxarray
+from shapely.geometry import mapping
+
+# Declare directories
+main_directory='.'
+
+# Define Plotting Info/Bounds
+dpi=300
+plt.style.use('dark_background')
+minLat = 33.75    
+maxLat = 36.75   
+minLon = -85 
+maxLon = -75
+
+####################
+# BEGIN PROGRAM
+####################
+start=time.time()
+
+# Read in Arguments 
+if len(sys.argv) < 6:
+    sys.exit("USAGE: <YYYY> <MM> <DD> <HH> <ELEMENT>\nExample: ipython plot_hrrr_nc-gridded.py 2021 08 25 12 T/Td")
+
+# Get Date Info
+year= "%04i" % int(sys.argv[1])
+month= "%02i" % int(sys.argv[2])
+day= "%02i" % int(sys.argv[3])
+hour = "%02i" % int(sys.argv[4])
+element= sys.argv[5]
+date = year+month+day
+init_name=year+month+day+': '+hour+'UTC'
+init_time=year+month+day+hour
+print("INIT: ",init_name)
+
+# Get Element Info
+if element == 'T':
+    var = 'TMP'
+    level = '2m_above_ground'
+    unit='$^\circ$F'
+    color_map='YlOrRd' 
+    outElement='2m Temperature'
+    bounds=np.arange(40,100,5)
+elif element =='Td':
+    var = 'DPT'
+    level = '2m_above_ground'
+    unit='$^\circ$F'
+    color_map='BrBG' 
+    outElement='2m Dew Point'
+    bounds=np.arange(40,100,5)
+else:
+    sys.exit("Only Available For Elements T/Td")
+
+#################################################
+# READ IN SHAPEFILE
+input_shapefile=main_directory+'/input_shapefile/NCSteMerc.shp'
+
+print("READ IN SHAPEFILE: ",input_shapefile)
+geo_shapefile = gpd.read_file(input_shapefile)
+
+# Get Projection
+projection=geo_shapefile.crs
+
+####################
+# PART ONE: Analysis
+#  i.e. T=0
+####################
+# Access HRRR from AWS ... projection dimensions are in url2
+# Note That 'anl' stands for Analysis (Initialization)
+fs = s3fs.S3FileSystem(anon=True)
+url1 = 's3://hrrrzarr/sfc/' + date + '/' + date + '_' + hour + 'z_anl.zarr/' + level + '/' + var + '/' + level
+url2 = 's3://hrrrzarr/sfc/' + date + '/' + date + '_' + hour + 'z_anl.zarr/' + level + '/' + var
+
+file1 = s3fs.S3Map(url1, s3=fs)
+file2 = s3fs.S3Map(url2, s3=fs)
+
+# Open Datasets
+print("\nOPEN ANALYSIS DATASETS")
+print ("\t",url1)
+print ("\t",url2)
+ds = xr.open_mfdataset([file1,file2], engine='zarr')
+
+# Define HRRR projection
+lon1 = -97.5
+lat1 = 38.5
+slat = 38.5
+projData= ccrs.LambertConformal(central_longitude=lon1,central_latitude=lat1,standard_parallels=[slat])
+
+# Get Data
+airTemp = ds[var].metpy.convert_units('degF')
+
+# Clip By Shapefile
+airTemp=airTemp.rio.set_spatial_dims(x_dim="projection_x_coordinate", y_dim="projection_y_coordinate", inplace=True)
+airTemp=airTemp.rio.write_crs(projData.proj4_params, inplace=True)
+#airTemp=airTemp.rio.reproject(projection)
+airTempC = airTemp.rio.clip(geo_shapefile.geometry.apply(mapping), geo_shapefile.crs, drop=False)
+
+# PLOT
+print("\nPLOT ANALYSIS")
+
+# Get Plotting Coordinates
+x = airTempC.projection_x_coordinate
+y = airTempC.projection_y_coordinate
+
+# Set Up ColorBar
+vmin=np.min(bounds)
+vmax=np.max(bounds)
+extend='both'
+norm = mcolors.BoundaryNorm(boundaries=bounds, ncolors=256)
+cm = plt.cm.ScalarMappable(cmap=color_map)
+cm.set_array(airTemp)
+cm.set_clim(vmin, vmax)
+
+# Set Up Figure
+# When trying to set projection, ccrs.Mercator() does not work
+outProjection=projData
+#outProjection=ccrs.Mercator()
+fig= plt.figure(num=1, figsize=(8,4), dpi=dpi, facecolor='w', edgecolor='k')
+ax = fig.add_axes([0, 0, 1, 1], projection=outProjection)
+ax.set_facecolor('black')
+ax.set_extent([minLon, maxLon, minLat, maxLat], crs=ccrs.PlateCarree())  
+
+# Add Boundaries
+#ax.add_feature(cfeature.COASTLINE,linewidth=0.5,edgecolor='white',facecolor='black',zorder=5)
+#ax.add_feature(cfeature.OCEAN,linewidth=0.5,edgecolor='white',facecolor='black',zorder=5)
+#ax.add_feature(cfeature.STATES,linewidth=0.5,edgecolor='white',facecolor='black',zorder=5)
+
+# Plot Data
+# The transform options do not work
+ax.pcolormesh(x, y, airTempC,cmap=color_map,norm=norm,vmin=vmin,vmax=vmax,zorder=8)
+#ax.pcolormesh(x, y, airTempC,transform=ccrs.PlateCarree(),cmap=color_map,norm=norm,vmin=vmin,vmax=vmax,zorder=8)
+#ax.pcolormesh(x, y, airTempC,transform=ccrs.Mercator(),cmap=color_map,norm=norm,vmin=vmin,vmax=vmax,zorder=8)
+
+# Add Colormap
+cax = fig.add_axes([0.1, -0.035, 0.8, 0.03])
+cbar=plt.colorbar(cm, cax=cax,boundaries=bounds,orientation='horizontal',extend=extend,spacing='uniform')
+cbar.ax.tick_params(labelsize=15)
+cbar.set_label(unit,size=15)
+
+# Add Titles
+plt.suptitle(outElement+' Analysis For '+init_name,fontsize=17,y=1.10)
+plt.annotate('Source: HRRR, init= '+init_name+'\nMade By Jared Rennie @jjrennie',xy=(1.045, -3.91), xycoords='axes fraction', fontsize=5,backgroundcolor='black',color='white',horizontalalignment='right', verticalalignment='bottom')
+
+# Save Figure
+plt.savefig(main_directory+"/hrrr-nc-grd_"+var.lower()+"_init-"+init_time+"-anls-00.png",bbox_inches='tight') 
+plt.clf()
+plt.close()
+
+####################
+# DONE
+####################
+print("DONE!")
+end=time.time()
+print ("Runtime: %8.1f seconds." % (end-start))
+sys.exit()
\ No newline at end of file