read_goes.py

   1 from netCDF4 import Dataset
   2 from datetime import datetime, timedelta
   3 import numpy as np
   4 import matplotlib.pyplot as plt
   5 from mpl_toolkits.basemap import Basemap
   6 from pyproj import Proj
   7 import sys
   8
   9 file=sys.argv[1]
  10 ncf=Dataset(file,'r')
  11
  12 # Seconds since 2000-01-01 12:00:00
  13 add_seconds=ncf.variables['t'][0]
  14 # Datetime of image scan
  15 dt=datetime(2000, 1, 1, 12)+timedelta(seconds=float(add_seconds.data))
  16 print 'dt = '
  17 print dt
  18
  19 # Detections in RGB array
  20 # Load the RGB arrays
  21 R=ncf.variables['CMI_C07'][:].data
  22 G=ncf.variables['CMI_C06'][:].data
  23 B=ncf.variables['CMI_C05'][:].data
  24 # Turn empty values into nans
  25 R[R==-1]=np.nan
  26 G[G==-1]=np.nan
  27 B[B==-1]=np.nan
  28 # Apply range limits for each channel (mostly important for Red channel)
  29 R=np.maximum(R, 273)
  30 R=np.minimum(R, 333)
  31 G=np.maximum(G, 0)
  32 G=np.minimum(G, 1)
  33 B=np.maximum(B, 0)
  34 B=np.minimum(B, .75)
  35 # Normalize each channel by the appropriate range of values (again, mostly important for Red channel)
  36 R=(R-273)/(333-273)
  37 G=(G-0)/(1-0)
  38 B=(B-0)/(.75-0)
  39 # Apply the gamma correction to Red channel.
  40 #   I was told gamma=0.4, but I get the right answer with gamma=2.5 (the reciprocal of 0.4)
  41 R=np.power(R, 2.5)
  42 # The final RGB array :)
  43 RGB=np.dstack([R, G, B])
  44 print 'RGB = '
  45 print RGB
  46
  47 # Geolocation
  48 # Satellite height
  49 sat_h=ncf.variables['goes_imager_projection'].perspective_point_height
  50 print 'sat_h = '
  51 print sat_h
  52 # Satellite longitude
  53 sat_lon=ncf.variables['goes_imager_projection'].longitude_of_projection_origin
  54 print 'sat_lon = '
  55 print sat_lon
  56 # Satellite sweep
  57 sat_sweep=ncf.variables['goes_imager_projection'].sweep_angle_axis
  58 print 'sat_sweep = '
  59 print sat_sweep
  60 # The projection x and y coordinates equals
  61 # the scanning angle (in radians) multiplied by the satellite height (http://proj4.org/projections/geos.html)
  62 X=ncf.variables['x'][:] * sat_h
  63 Y=ncf.variables['y'][:] * sat_h
  64 print 'X ='
  65 print X
  66 print 'Y ='
  67 print Y
  68
  69 # Visualize the way Brian Blaylock does
  70 p=Proj(proj='geos', h=sat_h, lon_0=sat_lon, sweep=sat_sweep)
  71 # Convert map points to latitude and longitude with the magic provided by Pyproj
  72 XX, YY=np.meshgrid(X, Y)
  73 lons, lats=p(XX, YY, inverse=True)
  74 lats[np.isnan(R)]=np.nan
  75 lons[np.isnan(R)]=np.nan
  76 # Make a new map object for the HRRR model domain map projection
  77 mH = Basemap(resolution='l', projection='lcc', area_thresh=5000, \
  78              width=1800*3000, height=1060*3000, \
  79              lat_1=38.5, lat_2=38.5, \
  80              lat_0=38.5, lon_0=-97.5)
  81 # Create a color tuple for pcolormesh
  82 rgb = RGB[:,:-1,:] # Using one less column is very important, else your image will be scrambled! (This is the stange nature of pcolormesh)
  83 colorTuple = rgb.reshape((rgb.shape[0] * rgb.shape[1]), 3) # flatten array, becuase that's what pcolormesh wants.
  84 colorTuple = np.insert(colorTuple, 3, 1.0, axis=1) # adding an alpha channel will plot faster, according to stackoverflow. Not sure why.
  85 # Now we can plot the GOES data on the HRRR map domain and projection
  86 plt.figure(figsize=[10,8])
  87 # The values of R are ignored becuase we plot the color in colorTuple, but pcolormesh still needs its shape.
  88 newmap = mH.pcolormesh(lons, lats, R, color=colorTuple, linewidth=0, latlon=True)
  89 newmap.set_array(None) # without this line the linewidth is set to zero, but the RGB colorTuple is ignored. I don't know why.
  90 mH.drawcoastlines(color='w')
  91 mH.drawcountries(color='w')
  92 mH.drawstates(color='w')
  93 plt.title('GOES-16 Fire Temperature', fontweight='semibold', fontsize=15)
  94 plt.title('%s' % dt.strftime('%H:%M UTC %d %B %Y'), loc='right')
  95 plt.show()
  96
  97 # location of South Sugarloaf fire
  98 l = {'latitude': 41.812,
  99      'longitude': -116.324}
 100 # Baltimore harbour
 101 #l = {'latitude': 39.2827513,
 102 #     'longitude':-76.6214623}
 103 # Buzzard fire
 104 #l = {'latitude': 33.724,
 105 #     'longitude': -108.538}
 106
 107 mZ = Basemap(projection='lcc', resolution='i',
 108             width=3E5, height=3E5,
 109             lat_0=l['latitude'], lon_0=l['longitude'],)
 110
 111 # Now we can plot the GOES data on a zoomed in map centered on the Sugarloaf wildfire
 112 plt.figure(figsize=[10, 10])
 113 newmap = mZ.pcolormesh(lons, lats, R, color=colorTuple, linewidth=0, latlon=True)
 114 newmap.set_array(None)
 115
 116 #mZ.scatter(l['longitude'], l['latitude'], marker='o', s=1000, facecolor='none', edgecolor='greenyellow') # circle fire area
 117 mZ.drawcoastlines(color='w')
 118 mZ.drawcountries(color='w')
 119 mZ.drawstates(color='w')
 120
 121 plt.title('GOES-16 Fire Temperature\n', fontweight='semibold', fontsize=15)
 122 plt.title('%s' % dt.strftime('%H:%M UTC %d %B %Y'), loc='right')
 123 plt.title('South Sugarloaf Fire', loc='left')
 124 plt.show()