Use Dave Tobin's AIRS channel properties file

airs2modis.py

@@ -3,8 +3,6 @@ import numpy
 import sys
 import time
 
-import matplotlib.pyplot as plt
-
 from modis_bright import modis_bright_from_airs
 from pyhdf.SD import SD, SDC
 
@@ -18,23 +16,13 @@ output_filename = sys.argv[2]
 # open our "helper" data file (see airs2modis_pre.py)
 helper_sd = SD('airs2modis.hdf')
 
-# from the given AIRS filename, determine its date and figure out
-# what generation of "channel properties" applies to it. this means
-# determining the most recent AIRS channel properties file that is
-# dated earlier or equal to the given granule
-st = time.strptime(airs_filename[:15], 'AIRS.%Y.%m.%d')
-airs_date = calendar.timegm(st)
-dates = helper_sd.select('Channel_Props_Date').get()
-date_idx = int((dates <= airs_date).nonzero()[0][-1])
-
 # read in the mask of good AIRS channels and apply it to the AIRS
 # radiances from our input granule. note that the mask is stored as
 # 8-bit integers and we need to change its data type to bool to
 # make numpy treat it as a mask and not an array of indices!
-airs_mask = helper_sd.select('Channel_Mask')[date_idx]
+airs_mask = helper_sd.select('Channel_Mask').get()
 airs_mask.dtype = numpy.bool8
 in_rads = SD(airs_filename).select('radiances').get()[:,:,airs_mask]
-num_good_chans = in_rads.shape[2]
 
 # bail if there are any remaining fill values in the so-called good
 # channels
@@ -48,11 +36,12 @@ out_rads = numpy.empty((16, 135, 90))
 # the weights from our helper file
 response_sds = helper_sd.select('Channel_Weights')
 for band_idx in range(16):
-    response = response_sds[date_idx,band_idx,:num_good_chans]
+    response = response_sds[band_idx,:]
     out_rads[band_idx] = (in_rads * response).sum(axis=2)
 
 # convert radiances to brightness temperatures
 bts = modis_bright_from_airs(out_rads.reshape((16, 135 * 90)))
+
 # write results out to HDF
 output_sd = SD(output_filename, SDC.WRITE | SDC.CREATE | SDC.TRUNC)
 bt_sds = output_sd.create('AIRS_Brightness_Temps', SDC.FLOAT32,

airs2modis_pre.py

-# This script uses a set of AIRS channel properties files and a MODIS
+# This script uses an AIRS channel properties files and a MODIS
 # spectral response file to produce an intermediate airs2modis.hdf
 # file. This intermediate file is then used as input by the
 # airs2modis module. The contents of the intermediate file include:
-#   - the times (in UNIX time) at which each channel properties file
-#     went into effect (Channel_Props_Date)
-#   - AIRS channel masks corresponding to the channels not flagged as
-#     bad, one for each channel properties file (Channel_Mask)
+#   - an AIRS channel mask corresponding to the channels not flagged
+#     as bad in the channel properties file (Channel_Mask)
 #   - channel weights that can be applied to the good channels from
-#     each mask, as interpolated from the MODIS SRF file
-#     (Channel_Weights)
+#     the mask, as interpolated from the MODIS SRF file for each
+#     MODIS band (Channel_Weights)
 
-import calendar
 import itertools
 import numpy
-import os
-import time
+import sys
 
 from pycdf import CDF
 from pyhdf.SD import SD, SDC
@@ -65,47 +61,38 @@ def read_airs_chan_data(filename):
 
     return mask, freqs[mask]
 
+# verify command line arguments
+if len(sys.argv) != 3:
+    print ('usage: python %s <airs_chan_props_file> <modis_srf_file>' %
+           sys.argv[0])
+    sys.exit(1)
+chan_props_file = sys.argv[1]
+srf_file = sys.argv[2]
+
 # read in the MODIS SRF data
-srf_data = read_modis_srf_data('modis_aqua.srf.nc')
+srf_data = read_modis_srf_data(srf_file)
 
 # the MODIS bands we are interested in getting SRF data for
 modis_bands = [20, 21, 22, 23, 24, 25, 27, 28,
                29, 30, 31, 32, 33, 34, 35, 36]
 num_bands = len(modis_bands)
 
-# get the list of AIRS channel properties files we'll be ingesting
-dirname = 'channel_properties_files'
-chan_prop_files = sorted(os.listdir(dirname))
-num_files = len(chan_prop_files)
-
-# initialize our HDF output
-sd = SD('airs2modis.hdf', SDC.WRITE | SDC.CREATE | SDC.TRUNC)
-date_sds = sd.create('Channel_Props_Date', SDC.INT32, (num_files,))
-mask_sds = sd.create('Channel_Mask', SDC.INT8, (num_files, 2378))
-response_sds = sd.create('Channel_Weights',
-                         SDC.FLOAT32,
-                         (num_files, num_bands, 2378))
-
-# loop through the channel propery files
-for file_num, filename in enumerate(chan_prop_files):
-
-    # determine the begin date for the file from its file
-    st = time.strptime(filename, 'L2.chan_prop.%Y.%m.%d.v9.5.1.txt')
-    date_sds[file_num] = calendar.timegm(st)
+# parse the channel props file for the mask and channel frequencies
+mask, freqs = read_airs_chan_data(chan_props_file)
 
-    # parse the file for the mask and channel frequencies
-    mask, freqs = read_airs_chan_data('%s/%s' % (dirname, filename))
-    mask_sds[file_num] = mask
+# determine how to weight the masked AIRS channels for each MODIS
+# band by interpolating the MODIS SRF data
+response = numpy.empty((num_bands, freqs.size), numpy.float64)
+for band_idx, band in enumerate(modis_bands):
+        response[band_idx,:] = numpy.interp(freqs,
+                                            srf_data[band].wavenumbers,
+                                            srf_data[band].weights,
+                                            0, 0)
+        response[band_idx,:] /= response[band_idx,:].sum()
 
-    # determine how to weight the masked AIRS channels for each MODIS
-    # band by interpolating the MODIS SRF data
-    num_good = freqs.size
-    for band_idx, band in enumerate(modis_bands):
-        response = numpy.empty((2378,), numpy.float32)
-        tmp = numpy.interp(freqs,
-                           srf_data[band].wavenumbers,
-                           srf_data[band].weights,
-                           0, 0)
-        response[:num_good] = tmp / tmp.sum()
-        response[num_good:] = 0.0
-        response_sds[file_num,band_idx] = response.reshape(1,2378)
+# write out the mask and weights to HDF
+sd = SD('airs2modis.hdf', SDC.WRITE | SDC.CREATE | SDC.TRUNC)
+mask_sds = sd.create('Channel_Mask', SDC.INT8, mask.shape)
+response_sds = sd.create('Channel_Weights', SDC.FLOAT64, response.shape)
+mask_sds[:] = mask
+response_sds[:] = response