In [1]:
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import datetime as dt
import cartopy
import cartopy.crs as ccrs
import Nio
%matplotlib inline
In [2]:
file='/scratch/scholar/m/mebaldwi/EAPS591/mrms/mrms_20170430-070000.grib2'
#open file(s)
gr = Nio.open_file(file,'r',format='grb2')
print gr
names =  gr.variables.keys()
print "Variable names:", names

#dimensions
dims = gr.dimensions
print "dimensions:",dims

#global attributes
attr = gr.attributes.keys()
print "global attributes",attr
Nio file:	mrms_20170430-070000.grib2.grb2
   global attributes:
   dimensions:
      ygrid_0 = 1059
      xgrid_0 = 1799
   variables:
      float SeamlessHSR_P0_L102_GLC0 [ ygrid_0, xgrid_0 ]
         center :	US NOAA Office of Oceanic and Atmospheric Research
         production_status :	Research products
         long_name :	Seamless Hybrid Scan Reflectivity with VPR Correction
         units :	dBZ
         _FillValue :	9.999e+20
         coordinates :	gridlat_0 gridlon_0
         grid_type :	Lambert Conformal can be secant or tangent, conical or bipolar
         parameter_discipline_and_category :	Multi-Radar/Multi-Sensor (MRMS) products, Category 8
         parameter_template_discipline_category_number :	[0, 209, 8, 8]
         level_type :	Specific altitude above mean sea level (m)
         level :	0
         forecast_time :	0
         forecast_time_units :	minutes
         initial_time :	04/30/2017 (07:00)
      float gridrot_0 [ ygrid_0, xgrid_0 ]
         long_name :	vector rotation angle
         GridType :	Lambert Conformal (secant, tangent, conical or bipolar)
         units :	radians
         formula_u :	Ugrid = cos(rot)*Uearth - sin(rot)*Vearth
         formula_v :	Vgrid = sin(rot)*Uearth + cos(rot)*Vearth
         note1 :	u and v components of vector quantities are resolved relative to earth
         note2 :	apply formulas to derive u and v components relative to grid
      float gridlat_0 [ ygrid_0, xgrid_0 ]
         corners :	[21.13812, 21.14055, 47.84219, 47.83862]
         long_name :	latitude
         grid_type :	Lambert Conformal (secant, tangent, conical or bipolar)
         units :	degrees_north
         Latin2 :	38.5
         Latin1 :	38.5
         Dy :	3
         Dx :	3
         Lov :	262.5
         Lo1 :	237.2805
         La1 :	21.13812
      float gridlon_0 [ ygrid_0, xgrid_0 ]
         corners :	[-122.7195, -72.28972, -60.91719, -134.0955]
         long_name :	longitude
         grid_type :	Lambert Conformal (secant, tangent, conical or bipolar)
         units :	degrees_east
         Latin2 :	38.5
         Latin1 :	38.5
         Dy :	3
         Dx :	3
         Lov :	262.5
         Lo1 :	237.2805
         La1 :	21.13812

Variable names: ['SeamlessHSR_P0_L102_GLC0', 'gridlat_0', 'gridrot_0', 'gridlon_0']
dimensions: {'xgrid_0': 1799, 'ygrid_0': 1059}
global attributes []
In [3]:
data = gr.variables['gridlat_0']
attrib = data.attributes.keys()
print attrib
dx = getattr(data,'Dx')
dy = getattr(data,'Dy')
print dx, dy
['Latin1', 'Latin2', 'corners', 'long_name', 'Lo1', 'Lov', 'Dx', 'Dy', 'units', 'La1', 'grid_type']
[ 3.] [ 3.]
In [7]:
#read in variables
obs = gr.variables["SeamlessHSR_P0_L102_GLC0"][:]
lats = gr.variables["gridlat_0"][:]
lons = gr.variables["gridlon_0"][:]
gr.close()
In [8]:
np.max(obs)
Out[8]:
57.5
In [9]:
np.min(obs)
Out[9]:
-999.0
In [10]:
# plot radar
fig = plt.figure()
fig.set_size_inches(12,8)
ax=plt.axes([0.,0.,1.,1.],projection=ccrs.PlateCarree())
ax.set_extent([-125.0,-66.5,22.5,52.5])
intervals = [0,5,10,15,20,25,30,35,40,45,50,55,60,65]
obsobj=plt.contourf(lons,lats,obs,intervals,cmap='jet',transform=ccrs.PlateCarree())
states_provinces = cartopy.feature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(states_provinces,edgecolor='gray')
colbar = plt.colorbar(obsobj,orientation='horizontal',ticks=intervals)
colbar.set_label('dBZ',rotation='horizontal',fontsize='medium')
titlestr = 'MRMS radar valid 20170430 0700 UTC '
plt.title(titlestr, fontsize='medium')
Out[10]:
<matplotlib.text.Text at 0x2b30bcfe0f50>
In [11]:
# plot radar missing area
fig = plt.figure()
fig.set_size_inches(12,8)
ax=plt.axes([0.,0.,1.,1.],projection=ccrs.PlateCarree())
ax.set_extent([-125.0,-66.5,22.5,52.5])
intervals = [-1000,-400,-300,-100,-50,-10,0]
obsobj=plt.contourf(lons,lats,obs,intervals,cmap='jet',transform=ccrs.PlateCarree())
states_provinces = cartopy.feature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(states_provinces,edgecolor='gray')
colbar = plt.colorbar(obsobj,orientation='horizontal',ticks=intervals)
colbar.set_label('dBZ',rotation='horizontal',fontsize='medium')
titlestr = 'MRMS radar valid 20170430 0700 UTC '
plt.title(titlestr, fontsize='medium')
Out[11]:
<matplotlib.text.Text at 0x2b30c30907d0>
In [12]:
hrrrfile='/scratch/scholar/m/mebaldwi/EAPS591/hrrr/hrrr.t06z.wrfsfcf01'
#open file(s)
gr2 = Nio.open_file(hrrrfile,'r',format='grb2')
print gr2
names =  gr2.variables.keys()
print "Variable names:", names

#dimensions
dims = gr2.dimensions
print "dimensions:",dims

#global attributes
attr = gr2.attributes.keys()
print "global attributes",attr
Nio file:	hrrr.t06z.wrfsfcf01.grb2
   global attributes:
   dimensions:
      ygrid_0 = 1059
      xgrid_0 = 1799
   variables:
      float REFD_P0_L103_GLC0 [ ygrid_0, xgrid_0 ]
         center :	US National Weather Service - NCEP (WMC)
         production_status :	Operational products
         long_name :	Reflectivity
         units :	dB
         _FillValue :	1e+20
         coordinates :	gridlat_0 gridlon_0
         grid_type :	Lambert Conformal can be secant or tangent, conical or bipolar
         parameter_discipline_and_category :	Meteorological products, Forecast radar imagery
         parameter_template_discipline_category_number :	[0, 0, 16, 195]
         level_type :	Specified height level above ground (m)
         level :	1000
         forecast_time :	1
         forecast_time_units :	hours
         initial_time :	04/30/2017 (06:00)
      float gridrot_0 [ ygrid_0, xgrid_0 ]
         long_name :	vector rotation angle
         GridType :	Lambert Conformal (secant, tangent, conical or bipolar)
         units :	radians
         formula_u :	Uearth = sin(rot)*Vgrid + cos(rot)*Ugrid
         formula_v :	Vearth = cos(rot)*Vgrid - sin(rot)*Ugrid
         note1 :	u and v components of vector quantities are resolved relative to grid
         note2 :	apply formulas to derive u and v components relative to earth
      float gridlat_0 [ ygrid_0, xgrid_0 ]
         corners :	[21.13812, 21.14055, 47.84219, 47.83862]
         long_name :	latitude
         grid_type :	Lambert Conformal (secant, tangent, conical or bipolar)
         units :	degrees_north
         Latin2 :	38.5
         Latin1 :	38.5
         Dy :	3
         Dx :	3
         Lov :	262.5
         Lo1 :	237.2805
         La1 :	21.13812
      float gridlon_0 [ ygrid_0, xgrid_0 ]
         corners :	[-122.7195, -72.28972, -60.91719, -134.0955]
         long_name :	longitude
         grid_type :	Lambert Conformal (secant, tangent, conical or bipolar)
         units :	degrees_east
         Latin2 :	38.5
         Latin1 :	38.5
         Dy :	3
         Dx :	3
         Lov :	262.5
         Lo1 :	237.2805
         La1 :	21.13812

Variable names: ['REFD_P0_L103_GLC0', 'gridlat_0', 'gridrot_0', 'gridlon_0']
dimensions: {'xgrid_0': 1799, 'ygrid_0': 1059}
global attributes []
In [13]:
#read in variables
fcst = gr2.variables["REFD_P0_L103_GLC0"][:]
hrrrlats = gr2.variables["gridlat_0"][:]
hrrrlons = gr2.variables["gridlon_0"][:]
In [14]:
gr2.close()
In [15]:
np.max(fcst)
Out[15]:
62.875
In [16]:
np.min(fcst)
Out[16]:
-10.0
In [17]:
# plot forecast
fig = plt.figure()
fig.set_size_inches(12,8)
ax=plt.axes([0.,0.,1.,1.],projection=ccrs.PlateCarree())
ax.set_extent([-125.0,-66.5,22.5,52.5])
intervals = [0,5,10,15,20,25,30,35,40,45,50,55,60,65]
fcstobj=plt.contourf(lons,lats,fcst,intervals,cmap='jet',transform=ccrs.PlateCarree())
states_provinces = cartopy.feature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(states_provinces,edgecolor='gray')
colbar = plt.colorbar(fcstobj,orientation='horizontal',ticks=intervals)
colbar.set_label('dBZ',rotation='horizontal',fontsize='medium')
titlestr = 'HRRR fcst radar valid 20170430 0700 UTC '
plt.title(titlestr, fontsize='medium')
Out[17]:
<matplotlib.text.Text at 0x2b30cc26ed10>
In [18]:
start_date = dt.datetime.strptime("2017-04-30_06:00:00","%Y-%m-%d_%H:%M:%S")
duration=18
In [19]:
print np.mean(obs[obs>=-10])
print np.mean(fcst[obs>=-10])
18.9446346818
10.8223
In [20]:
[imax,jmax]=fcst.shape
thresh=20.
asum=np.zeros((imax,jmax))
bsum=np.zeros((imax,jmax))
csum=np.zeros((imax,jmax))
dsum=np.zeros((imax,jmax))
print imax,jmax
for i in np.arange(duration):
    fhr=i+1
    forecast_hour='%02d'%int(fhr)
    valid_time=start_date + dt.timedelta(hours=fhr)
    valid_date=valid_time.strftime("%Y%m%d-%H%M%S")
    #print valid_time
    hrrrfile='/scratch/scholar/m/mebaldwi/EAPS591/hrrr/hrrr.t06z.wrfsfcf'+forecast_hour
    mrmsfile='/scratch/scholar/m/mebaldwi/EAPS591/mrms/mrms_' + valid_date + '.grib2'
    #print hrrrfile
    #print mrmsfile
    gr = Nio.open_file(mrmsfile,'r',format='grb2')
    obsgrid = gr.variables["SeamlessHSR_P0_L102_GLC0"][:]
    gr.close()
    gr2 = Nio.open_file(hrrrfile,'r',format='grb2')
    fcstgrid = gr2.variables["REFD_P0_L103_GLC0"][:]
    gr2.close()
    obs=np.ma.masked_where(obsgrid<-99.,obsgrid)
    fcst=np.ma.masked_where(obsgrid<-99.,fcstgrid)
    asum=asum+np.logical_and(fcst>thresh,obs>thresh)
    bsum=bsum+np.logical_and(fcst>thresh,obs<=thresh)
    csum=csum+np.logical_and(fcst<=thresh,obs>thresh)
    dsum=dsum+np.logical_and(fcst<=thresh,obs<=thresh)
    
a=np.sum(asum)
b=np.sum(bsum)
c=np.sum(csum)
d=np.sum(dsum)
print a,b,c,d
1059 1799
566621.0 996554.0 484539.0 25829894.0
In [21]:
# plot forecast
fig = plt.figure()
fig.set_size_inches(12,8)
ax=plt.axes([0.,0.,1.,1.],projection=ccrs.PlateCarree())
ax.set_extent([-125.0,-66.5,22.5,52.5])
aobj=plt.contourf(lons,lats,asum,cmap='jet',transform=ccrs.PlateCarree())
states_provinces = cartopy.feature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(states_provinces,edgecolor='gray')
colbar = plt.colorbar(aobj,orientation='horizontal')
colbar.set_label('asum',rotation='horizontal',fontsize='medium')
titlestr = 'HRRR vs MRMS 18h fcst 2017043007-2017050100 '
plt.title(titlestr, fontsize='medium')
Out[21]:
<matplotlib.text.Text at 0x2b30ccaad410>
In [22]:
print (a/(a+b+c))
0.276709052143
In [ ]: