read_goes.py

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