Added modis2airs that distinguishes mirror / detector

modis2airs_striping.py

+import numpy
+import pyhdf.VS
+import sys
+
+from flags import *
+from modis_bright import modis_bright
+from pyhdf.HDF import HDF
+from pyhdf.SD import SD, SDC
+
+# helper function for reading in MODIS radiances. returns a numpy
+# array with dimensions:
+#   - 16 (MODIS band, 20-25,27-36)
+#   - 2030 or 2040 (MODIS along track index)
+#   - 1354 (MODIS across track index)
+def read_modis_radiances(filename):
+
+    sd = SD(filename)
+    sds = sd.select('EV_1KM_Emissive')
+    scales = numpy.array(sds.radiance_scales).reshape((16, 1, 1))
+    offsets = numpy.array(sds.radiance_offsets).reshape((16, 1, 1))
+    raw = sds.get()
+    bad_idxs = (raw >= 32768)
+    rads = scales * (raw - offsets)
+    rads[bad_idxs] = -9999.0
+    return rads
+
+
+# helper function for reading in MODIS mirror sides. returns a list
+# that has a 0/1 value for each scan in the granule (so either 203
+# or 204 values)
+def read_modis_mirror_sides(filename):
+
+    hdf = HDF(filename)
+    vs = hdf.vstart()
+    vd = vs.attach('Level 1B Swath Metadata')
+    vd.setfields('Mirror Side')
+    return [x[0] for x in vd[:]]
+
+# verify command line arguments
+if len(sys.argv) != 4:
+    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]
+output_file = sys.argv[3]
+
+# read in the MODIS radiances
+in_rads = read_modis_radiances(modis_file)
+
+# read in the MODIS mirror sides
+mirror_sides = read_modis_mirror_sides(modis_file)
+
+# and the collocation data
+coll_sd = SD(coll_file)
+num_colls = coll_sd.select('Number_Of_MODIS_FOV').get()
+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, 2, 10), numpy.float32)
+rad_std_arr = numpy.empty((16, 135, 90, 2, 10), numpy.float32)
+bt_arr = numpy.empty((16, 135, 90, 2, 10), numpy.float32)
+std_bt_arr = numpy.empty((16, 135, 90, 2, 10), numpy.float32)
+n_arr = numpy.empty((16, 135, 90, 2, 10), numpy.int32)
+flags_arr = numpy.zeros((16, 135, 90, 2, 10), numpy.int8)
+
+# loop over each AIRS pixel
+for airs_row in range(135):
+    for airs_col in range(90):
+
+        # loop through the collocations for this AIRS pixel,
+        # separating them based on mirror and detector number
+        idx_lists = [[([], []) for j in range(10)] for i in range(2)]
+        for i in range(num_colls[airs_row,airs_col]):
+            along_idx = modis_along_idxs[airs_row,airs_col,i]
+            across_idx = modis_across_idxs[airs_row,airs_col,i]
+            mirror = mirror_sides[along_idx / 10]
+            detector = along_idx % 10
+            idx_lists[mirror][detector][0].append(along_idx)
+            idx_lists[mirror][detector][1].append(across_idx)
+
+        # now loop through the mirror / detector combinations,
+        # compiling the collocated radiances
+        for mirror in range(2):
+            for detector in range(10):
+
+                # avoid excessive indexing in the code below
+                rads_mean = rad_mean_arr[:,airs_row,airs_col,mirror,detector]
+                rads_std = rad_std_arr[:,airs_row,airs_col,mirror,detector]
+                bts = bt_arr[:,airs_row,airs_col,mirror,detector]
+                std_bts = std_bt_arr[:,airs_row,airs_col,mirror,detector]
+                n = n_arr[:,airs_row,airs_col,mirror,detector]
+                flags = flags_arr[:,airs_row,airs_col,mirror,detector]
+
+                # check for special case of zero collocations
+                num = len(idx_lists[mirror][detector][0])
+                if num == 0:
+                    rads_mean[:] = -9999.0
+                    rads_std[:] = -9999.0
+                    bts[:] = -9999.0
+                    std_bts[:] = -9999.0
+                    n[:] = 0
+                    flags |= MISSING
+                    continue
+
+                # pull out the collocated radiances
+                along_idxs = idx_lists[mirror][detector][0]
+                across_idxs = idx_lists[mirror][detector][1]
+                rads = in_rads[:,along_idxs,across_idxs]
+
+                # 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)
+
+                # 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)
+
+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_bt_sds.setfillvalue(-9999.0)
+
+rad_mean_sds[:] = rad_mean_arr
+rad_std_sds[:] = rad_std_arr
+bt_sds[:] = bt_arr
+std_bt_sds[:] = std_bt_arr
+n_sds[:] = n_arr
+flags_sds[:] = flags_arr