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='/depot/mebaldwi/data/archive/grib/mrms/20170430/mrms_20170430-180000.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-180000.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 (18: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 [4]:
#read in variables
obs = gr.variables["SeamlessHSR_P0_L102_GLC0"][:]
lats = gr.variables["gridlat_0"][:]
lons = gr.variables["gridlon_0"][:]
gr.close()
In [5]:
np.max(obs)
Out[5]:
54.0
In [6]:
np.min(obs)
Out[6]:
-999.0
In [7]:
# plot radar
fig = plt.figure()
fig.set_size_inches(12,8)
ax=plt.axes([0.,0.,1.,1.],projection=ccrs.PlateCarree())
ax.set_extent([-95.0,-82.0,28.0,39.0])
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 1800 UTC '
plt.title(titlestr, fontsize='medium')
Out[7]:
<matplotlib.text.Text at 0x2aec8541acd0>
In [8]:
# 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[8]:
<matplotlib.text.Text at 0x2aec8b56c950>
In [9]:
hrrrfile='/depot/mebaldwi/data/archive/grib/hrrr/2017043006/hrrr.t06z.wrfsfcf12'
#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.wrfsfcf12.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 :	12
         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 [10]:
#read in variables
fcst = gr2.variables["REFD_P0_L103_GLC0"][:]
hrrrlats = gr2.variables["gridlat_0"][:]
hrrrlons = gr2.variables["gridlon_0"][:]
gr2.close()
In [11]:
gr2.close()
In [12]:
np.max(fcst)
Out[12]:
57.625
In [13]:
np.min(fcst)
Out[13]:
-10.0
In [14]:
# 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])
ax.set_extent([-95.0,-82.0,28.0,39.0])
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 12h fcst REFD 1km AGL valid 20170430 1800 UTC '
plt.title(titlestr, fontsize='medium')
Out[14]:
<matplotlib.text.Text at 0x2aec9433c610>
In [20]:
start_date = dt.datetime.strptime("2017-03-25_00:00:00","%Y-%m-%d_%H:%M:%S")
duration=18
In [21]:
print np.mean(obs[obs>=-10])
print np.mean(fcst[obs>=-10])
16.6151567226
6.96322647443
In [ ]:
[imax,jmax]=fcst.shape
numcyc=888
duration=18
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 j in np.arange(numcyc):
    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")
        yymmdd_valid=valid_time.strftime("%Y%m%d")
        cycle_str=start_date.strftime("%H")
        yymmddhh_start=start_date.strftime("%Y%m%d%H")
        hrrrfile='/depot/mebaldwi/data/archive/grib/hrrr/'+yymmddhh_start+'/hrrr.t'+cycle_str+'z.wrfsfcf'+forecast_hour
        mrmsfile='/depot/mebaldwi/data/archive/grib/mrms/'+yymmdd_valid+'/mrms_' + valid_date + '.grib2'
        #file names with variable dates and times
        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 
In [18]:
print mrmsfile
/depot/mebaldwi/data/archive/grib/mrms/20170501/mrms_20170501-000000.grib2
In [19]:
# 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 fcst and obs > threshold'
plt.title(titlestr, fontsize='medium')
Out[19]:
<matplotlib.text.Text at 0x2aec7edd87d0>
In [22]:
print (a/(a+b+c))
0.276709052143
In [ ]: