Added validity checking of input data, among other changes

airs2modis.py

 import calendar
+import flags
 import numpy
 import sys
 import time
@@ -24,26 +25,47 @@ airs_mask = helper_sd.select('Channel_Mask').get()
 airs_mask.dtype = numpy.bool8
 in_rads = SD(airs_filename).select('radiances').get()[:,:,airs_mask]
 
-# bail if there are any remaining fill values in the so-called good
-# channels
-if (in_rads == -9999.0).any():
-    raise Exception('ERROR: radiances contain fill values')
-
-# set up array for output radiances
-out_rads = numpy.empty((16, 135, 90))
+# set up arrays for output radiances and flags
+out_rads = numpy.empty((16, 135, 90), numpy.float64)
+out_flags = numpy.zeros((16, 135, 90), numpy.int8)
 
 # loop over the MODIS bands, multiplying the AIRS spectra by the
 # the weights from our helper file
-response_sds = helper_sd.select('Channel_Weights')
+responses = helper_sd.select('Channel_Weights').get()
 for band_idx in range(16):
-    response = response_sds[band_idx,:]
-    out_rads[band_idx] = (in_rads * response).sum(axis=2)
+    out_rads[band_idx] = (in_rads * responses[band_idx,:]).sum(axis=2)
 
 # convert radiances to brightness temperatures
-bts = modis_bright_from_airs(out_rads.reshape((16, 135 * 90)))
+out_bts = modis_bright_from_airs(out_rads)
+
+# check for fill data that may have affected our results and replace
+# any such results with fill
+response_mask = (responses > 0.0)[:,numpy.newaxis,numpy.newaxis,:]
+invalid_mask = ((in_rads == -9999.0)[numpy.newaxis,:,:,:] &
+                response_mask).any(axis=3)
+out_rads[invalid_mask] = -9999.0
+out_bts[invalid_mask] = -9999.0
+out_flags[invalid_mask] |= flags.INVALID
+
+# check for negative values below what would be expected based on the
+# instrument's noise characteristics
+low_values = -1 * helper_sd.select('Channel_Noise_Rating').get()
+low_mask = ((in_rads < low_values)[numpy.newaxis,:,:,:] &
+            response_mask).any(axis=3)
+out_flags[low_mask] |= flags.LOW
+
+# check for values that seem too high (i.e. would yield a BT greater
+# than 350K)
+high_values = helper_sd.select('Channel_350K_Radiance').get()
+high_mask = ((in_rads > high_values)[numpy.newaxis,:,:,:] &
+             response_mask).any(axis=3)
+out_flags[high_mask] |= flags.HIGH
 
 # 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,
-                          (16, 135, 90))
-bt_sds[:] = numpy.array(bts, numpy.float32)
+rad_sds = output_sd.create('AIRS_Radiance', SDC.FLOAT64, out_rads.shape)
+bt_sds = output_sd.create('AIRS_Brightness_Temp', SDC.FLOAT64, out_bts.shape)
+flags_sds = output_sd.create('AIRS_Flags', SDC.INT8, out_flags.shape)
+rad_sds[:] = out_rads
+bt_sds[:] = out_bts
+flags_sds[:] = out_flags

airs2modis_pre.py

@@ -12,6 +12,7 @@ import itertools
 import numpy
 import sys
 
+from bright import inverse_bright_wavenumber
 from pyhdf.SD import SD, SDC
 
 # returns a dictionary of objects providing spectral response
@@ -36,10 +37,13 @@ def read_modis_srf_data(filename):
 
     return srf_data
 
-# returns a tuple with two pieces of information from the given AIRS
-# channel properties file:
-#   1. an array with the channel frequencies
-#   2. an array mask to filter out channels flagged as bad
+# returns a tuple with three pieces of information from the given
+# AIRS channel properties file:
+#   1. an array mask to filter out channels flagged as bad
+#   2. an array with the channel frequencies
+#   3. an array with the channel noise equivalent delta temperatures
+#      at 250K
+#   
 def read_airs_chan_data(filename):
 
     # loop through lines in the text file:
@@ -54,11 +58,13 @@ def read_airs_chan_data(filename):
     size = len(lines) - offset
     mask = numpy.empty((size,), numpy.bool)
     freqs = numpy.empty((size,), numpy.float64)
+    nedts = numpy.empty((size,), numpy.float64)
     for i in range(offset, len(lines)):
         mask[i - offset] = numpy.bool(lines[i][76] == '0')
         freqs[i - offset] = numpy.float64(lines[i][6:14])
+        nedts[i - offset] = numpy.float64(lines[i][26:32])
 
-    return mask, freqs[mask]
+    return mask, freqs[mask], nedts[mask]
 
 # verify command line arguments
 if len(sys.argv) != 3:
@@ -77,7 +83,7 @@ modis_bands = [20, 21, 22, 23, 24, 25, 27, 28,
 num_bands = len(modis_bands)
 
 # parse the channel props file for the mask and channel frequencies
-mask, freqs = read_airs_chan_data(chan_props_file)
+mask, freqs, nedts = read_airs_chan_data(chan_props_file)
 
 # determine how to weight the masked AIRS channels for each MODIS
 # band by interpolating the MODIS SRF data
@@ -89,9 +95,24 @@ for band_idx, band in enumerate(modis_bands):
                                             0, 0)
         response[band_idx,:] /= response[band_idx,:].sum()
 
+# convert the channel noise equivalent delta temperatures (NeDTs) to
+# a value in terms of radiance
+bt250 = 250.0 * numpy.ones((nedts.size,), numpy.float64)
+noise = (inverse_bright_wavenumber(freqs, bt250 + nedts) -
+         inverse_bright_wavenumber(freqs, bt250))
+
+# figure out what radiance would produce a BT of 350K for each
+# channel
+bt350 = 350.0 * numpy.ones((freqs.size,), numpy.float64)
+r350 = inverse_bright_wavenumber(freqs, bt350)
+
 # 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)
+noise_sds = sd.create('Channel_Noise_Rating', SDC.FLOAT64, noise.shape)
+r350_sds = sd.create('Channel_350K_Radiance', SDC.FLOAT64, r350.shape)
 mask_sds[:] = mask
 response_sds[:] = response
+noise_sds[:] = noise
+r350_sds[:] = r350

bright.py

+# General brightness temperature routines
+
+import numpy
+
+# Planck's constant (Joule seconds)
+h = 6.6260755e-34
+
+# speed of light in a vacuum (meters per second)
+c = 2.9979246e+8
+
+# Boltzmann constant (Joules per Kelvin)
+k = 1.380658e-23
+
+# derived constants
+c1 = 2.0 * h * c * c
+c2 = h * c / k
+
+# compute brightness temperature with inputs of wavelength in
+# microns and radiance in Watts per square meter per steradian
+# per micron
+def bright_wavelength(wavelength, radiance):
+
+    wavelength = wavelength * 1e-6
+    radiance = radiance * 1e6
+    return (c2 /
+            (wavelength * numpy.log(1.0 + c1 /
+                                          (radiance *
+                                           numpy.power(wavelength, 5)))))
+
+# compute brightness temperature with inputs of wavenumber in
+# inverse centimeters and radiance in milliWatts per square meter per
+# steradian per inverse centimeter
+def bright_wavenumber(wavenumber, radiance):
+
+    wavenumber = wavenumber * 1e2
+    radiance = radiance * 1e-5
+    return (c2 *
+            wavenumber /
+            numpy.log(1.0 + ((c1 * numpy.power(wavenumber, 3)) / radiance)))
+
+# compute radiance in milliWatts per meters squared per steradians
+# per inverse centimeters given inputs of wavenumber in inverse
+# centimeters and brightness temperature in Kelvin
+def inverse_bright_wavenumber(wavenumber, bt):
+
+    wavenumber = wavenumber * 1e2
+    rad = ((c1 * numpy.power(wavenumber, 3)) /
+           (numpy.exp(c2 * wavenumber / bt) - 1.0))
+    rad *= 1e5
+    return rad

flags.py

+# Module with values for flags output in the intercal system
+
+MISSING = 1
+INVALID = 2
+LOW = 4
+HIGH = 8

modis2airs.py

 import numpy
 import sys
 
+from flags import *
 from modis_bright import modis_bright
 from pyhdf.SD import SD, SDC
 
@@ -16,14 +17,15 @@ def read_modis_radiances(filename):
     scales = numpy.array(sds.radiance_scales).reshape((16, 1, 1))
     offsets = numpy.array(sds.radiance_offsets).reshape((16, 1, 1))
     raw = sds.get()
-    if (raw >= 32768).any():
-        raise Exception('ERROR: radiances contain invalid data')
+    bad_idxs = (raw >= 32768)
     rads = scales * (raw - offsets)
+    rads[bad_idxs] = -9999.0
     return rads
     
 # verify command line arguments
 if len(sys.argv) != 4:
-    print 'usage: %s <modis_file> <collocation_file> <output_file>' % sys.argv[0]
+    print ('usage: %s <modis_file> <collocation_file> <output_file>' %
+           sys.argv[0])
     sys.exit(1)
 modis_file = sys.argv[1]
 coll_file = sys.argv[2]
@@ -39,42 +41,88 @@ modis_along_idxs = coll_sd.select('MODIS_Along_Track_IDX').get()
 modis_across_idxs = coll_sd.select('MODIS_Across_Track_IDX').get()
 
 # set up numpy arrays for output
+rad_mean_arr = numpy.empty((16, 135, 90), numpy.float32)
+rad_std_arr = numpy.empty((16, 135, 90), numpy.float32)
 bt_arr = numpy.empty((16, 135, 90), numpy.float32)
-std_arr = numpy.empty((16, 135, 90), numpy.float32)
-n_arr = numpy.empty((135, 90), numpy.int32)
+std_bt_arr = numpy.empty((16, 135, 90), numpy.float32)
+n_arr = numpy.empty((16, 135, 90), numpy.int32)
+flags_arr = numpy.zeros((16, 135, 90), numpy.int8)
 
 # loop over each AIRS pixel
 for airs_row in range(135):
     for airs_col in range(90):
 
+        # avoid excessive indexing in the code below
+        rads_mean = rad_mean_arr[:,airs_row,airs_col]
+        rads_std = rad_std_arr[:,airs_row,airs_col]
+        bts = bt_arr[:,airs_row,airs_col]
+        std_bts = std_bt_arr[:,airs_row,airs_col]
+        n = n_arr[:,airs_row,airs_col]
+        flags = flags_arr[:,airs_row,airs_col]
+
         # use the collocation data to pull out overlapping MODIS
         # radiances (skip this AIRS FOV if there are none)
         num = num_colls[airs_row,airs_col]
         if num == 0:
-            bt_arr[:,airs_row,airs_col] = -9999.0
-            std_arr[:,airs_row,airs_col] = -9999.0
-            n_arr[airs_row,airs_col] = 0
+            rads_mean[:] = -9999.0
+            rads_std[:] = -9999.0
+            bts[:] = -9999.0
+            std_bts[:] = -9999.0
+            n[:] = 0
+            flags |= MISSING
             continue
         along_idxs = modis_along_idxs[airs_row,airs_col,:num]
         across_idxs = modis_across_idxs[airs_row,airs_col,:num]
         rads = in_rads[:,along_idxs,across_idxs]
 
-        # convert the collocated radiances to brightness temperature
-        bts = modis_bright(rads)
+        # if there's any invalid data amongst the collocated MODIS
+        # pixels, filter it out and set the INVALID flag for the
+        # affected bands
+        bad_idxs = (rads == -9999.0)
+        flags[bad_idxs.any(axis=1)] |= INVALID
+
+        # determine how many pixels are being used for each band
+        n[:] = num - bad_idxs.sum(axis=1)
+
+        # get mean radiances (bad measurements are set to zero so
+        # they have no effect on the result)
+        # FIXME: use masked arrays for this instead?
+        rads[bad_idxs] = 0.0
+        rads_mean[:] = rads.sum(axis=1) / n
+
+        # get standard deviation (bad measurements are cancelled out
+        # so they have no effect on the result)
+        diffs = rads - rads_mean.reshape((16, 1))
+        diffs[bad_idxs] = 0.0
+        rads_std[:] = numpy.sqrt((diffs * diffs).sum(axis=1)) / n
+
+        # now convert the radiances to brightness temperature
+        bts[:] = modis_bright(rads_mean)
 
-        # save output data: MODIS radiance mean and standard
-        # deviation, and number of MODIS pixels used
-        bt_arr[:,airs_row,airs_col] = bts.mean(axis=1)
-        std_arr[:,airs_row,airs_col] = bts.std(axis=1)
-        n_arr[airs_row,airs_col] = num
+        # to get a feel for the standard deviation in terms of
+        # brightness temperature, see what BT difference results from
+        # a radiance one std above the mean
+        std_bts[:] = modis_bright(rads_mean + rads_std) - bts
 
 # write our results to HDF
+
 sd = SD(output_file, SDC.WRITE | SDC.CREATE | SDC.TRUNC)
-bt_sds = sd.create('MODIS_Brightness_Temp', SDC.FLOAT32, (16, 135, 90))
-std_sds = sd.create('MODIS_Brightness_Temp_Std', SDC.FLOAT32, (16, 135, 90))
-n_sds = sd.create('MODIS_Brightness_Temp_N', SDC.INT32, (135, 90))
+
+rad_mean_sds = sd.create('MODIS_Radiance', SDC.FLOAT32, rad_mean_arr.shape)
+rad_std_sds = sd.create('MODIS_Radiance_Std', SDC.FLOAT32, rad_std_arr.shape)
+bt_sds = sd.create('MODIS_Brightness_Temp', SDC.FLOAT32, bt_arr.shape)
+std_bt_sds = sd.create('MODIS_Std_Brightness', SDC.FLOAT32, std_bt_arr.shape)
+n_sds = sd.create('MODIS_N', SDC.INT32, n_arr.shape)
+flags_sds = sd.create('MODIS_Flags', SDC.INT8, flags_arr.shape)
+
+rad_mean_sds.setfillvalue(-9999.0)
+rad_std_sds.setfillvalue(-9999.0)
 bt_sds.setfillvalue(-9999.0)
-std_sds.setfillvalue(-9999.0)
+std_bt_sds.setfillvalue(-9999.0)
+
+rad_mean_sds[:] = rad_mean_arr
+rad_std_sds[:] = rad_std_arr
 bt_sds[:] = bt_arr
-std_sds[:] = std_arr
+std_bt_sds[:] = std_bt_arr
 n_sds[:] = n_arr
+flags_sds[:] = flags_arr

modis_bright.py

@@ -3,26 +3,28 @@
 
 import numpy
 
-# centroid wavenumbers for bands 20-25,27-36 (inverse meters)
-centroid_wn = numpy.array([264740.89,
-                           251176.03,
-                           251790.84,
-                           246244.17,
-                           224829.57,
-                           220954.64,
-                           147426.20,
-                           136162.65,
-                           116962.56,
-                           102873.95,
-                           90768.13,
-                           83084.11,
-                           74829.77,
-                           73077.65,
-                           71820.94,
-                           70350.07]).reshape((16, 1))
-
-# centroid wavelengths (meters)
-centroid_wl = 1.0 / centroid_wn
+from bright import bright_wavelength, bright_wavenumber
+
+# centroid wavenumbers for bands 20-25,27-36 (inverse centimeters)
+centroid_wn = numpy.array([2647.4089,
+                           2511.7603,
+                           2517.9084,
+                           2462.4417,
+                           2248.2957,
+                           2209.5464,
+                           1474.2620,
+                           1361.6265,
+                           1169.6256,
+                           1028.7395,
+                           907.6813,
+                           830.8411,
+                           748.2977,
+                           730.7765,
+                           718.2094,
+                           703.5007]).reshape((16, 1))
+
+# centroid wavelengths (microns)
+centroid_wl = 1.0e4 / centroid_wn
 
 # temperature correction slopes for bands 20-25,27-36
 temp_corr_slope = numpy.array([0.99934,
@@ -60,37 +62,45 @@ temp_corr_int = numpy.array([0.48007,
                              0.01912,
                              0.01793]).reshape((16, 1))
 
-# Planck's constant (Joule seconds)
-h = 6.6260755e-34
-
-# speed of light in a vacuum (meters per second)
-c = 2.9979246e+8
-
-# Boltzmann constant (Joules per Kelvin)
-k = 1.380658e-23
+# helper function to make an array 2-dimensional. the first dimension
+# is always left the same size. if there is only 1 dimension, a 2nd
+# dimension of size 1 is added. if there are more than 2 dimensions,
+# the first dimension is left alone and the rest are smushed together
+def make_2d(arr):
 
-# derived constants
-c1 = 2.0 * h * c * c
-c2 = h * c / k
+    if arr.ndim == 1:
+        return arr.reshape((arr.size, 1))
+    else:
+        return arr.reshape((arr.shape[0],
+                            numpy.array(arr.shape[1:]).prod()))
 
 # convert MODIS radiances (in Watts per square meter per steradian
 # per micron) to brightness temperature
 def modis_bright(rad):
 
-    bt = (c2 /
-          (centroid_wl * numpy.log(1.0 + c1 /
-                                         (1.0e6 *
-                                          rad *
-                                          numpy.power(centroid_wl, 5)))))
-    return (bt - temp_corr_int) / temp_corr_slope
+    # we need to operate on a 16xN array
+    in_shape = rad.shape
+    rad = make_2d(rad)
+
+    # apply the Planck then the band-specific correction
+    bt = bright_wavelength(centroid_wl, rad)
+    bt = (bt - temp_corr_int) / temp_corr_slope
+
+    # return an array shaped like the one passed in
+    return bt.reshape(in_shape)
 
 # convert AIRS radiances (in milliWatts per square meter per
 # steradian per inverse centimeter) stored in terms of MODIS band
 # responses to brightness temperature
 def modis_bright_from_airs(rad):
 
-    bt = (c2 *
-          centroid_wn /
-          numpy.log(1.0 + ((c1 * numpy.power(centroid_wn, 3)) /
-                           (1.0e-5 * rad))))
-    return (bt - temp_corr_int) / temp_corr_slope
+    # we need to operate on a 16xN array
+    in_shape = rad.shape
+    rad = make_2d(rad)
+
+    # apply the Planck then the band-specific correction
+    bt = bright_wavenumber(centroid_wn, rad)
+    bt = (bt - temp_corr_int) / temp_corr_slope
+
+    # return an array shaped like the one passed in
+    return bt.reshape(in_shape)