Added modis2airs_collect.py

This allows us to consider everything in terms of AIRS
spatial resolution.

airs2modis.py

@@ -2,7 +2,6 @@ import calendar
 import flags
 import numpy
 import sys
-import time
 
 from modis_bright import modis_bright_from_airs
 from pyhdf.SD import SD, SDC

flags.py

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

modis2airs.py

@@ -2,7 +2,6 @@ import numpy
 import sys
 
 from flags import *
-from modis_bright import modis_bright
 from pyhdf.SD import SD, SDC
 
 # helper function for reading in MODIS radiances. returns a numpy
@@ -41,11 +40,9 @@ 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_bt_arr = numpy.empty((16, 135, 90), numpy.float32)
 n_arr = numpy.empty((16, 135, 90), numpy.int32)
+rad_sum_arr = numpy.empty((16, 135, 90), numpy.float64)
+rad_squared_sum_arr = numpy.empty((16, 135, 90), numpy.float64)
 flags_arr = numpy.zeros((16, 135, 90), numpy.int8)
 
 # loop over each AIRS pixel
@@ -53,23 +50,18 @@ 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]
+        rads_sum = rad_sum_arr[:,airs_row,airs_col]
+        rads_squared_sum = rad_squared_sum_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:
-            rads_mean[:] = -9999.0
-            rads_std[:] = -9999.0
-            bts[:] = -9999.0
-            std_bts[:] = -9999.0
             n[:] = 0
-            flags |= MISSING
+            rads_sum[:] = 0.0
+            rads_squared_sum[:] = 0.0
             continue
         along_idxs = modis_along_idxs[airs_row,airs_col,:num]
         across_idxs = modis_across_idxs[airs_row,airs_col,:num]
@@ -81,48 +73,28 @@ for airs_row in range(135):
         bad_idxs = (rads == -9999.0)
         flags[bad_idxs.any(axis=1)] |= INVALID
 
-        # determine how many pixels are being used for each band
+        # determine how many pixels are being used per 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)
+        # get radiance sum and squared sum (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
+        rads_sum[:] = rads.sum(axis=1)
+        rads_squared_sum[:] = (rads * rads).sum(axis=1)
 
-        # 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
+# write 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)
+rad_sum_sds = sd.create('MODIS_Radiance_Sum', SDC.FLOAT64, rad_sum_arr.shape)
+rad_squared_sum_sds = sd.create('MODIS_Radiance_Squared_Sum',
+                                SDC.FLOAT64,
+                                rad_squared_sum_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
+rad_sum_sds[:] = rad_sum_arr
+rad_squared_sum_sds[:] = rad_squared_sum_arr
 flags_sds[:] = flags_arr

modis2airs_collect.py

+import numpy
+import sys
+
+from modis_bright import modis_bright
+from pyhdf.SD import SD, SDC
+
+# verify command line arguments
+if len(sys.argv) < 2:
+    print 'usage: python %s <modis_file> [...] <output_file>' % sys.argv[0]
+    sys.exit(1)
+first_file = sys.argv[1]
+other_files = sys.argv[2:-1]
+output_file = sys.argv[-1]
+
+# read in data from the first intermediate file
+sd = SD(first_file)
+rads_sum = sd.select('MODIS_Radiance_Sum').get()
+rads_squared_sum = sd.select('MODIS_Radiance_Squared_Sum').get()
+n = sd.select('MODIS_N').get()
+flags = sd.select('MODIS_Flags').get()
+
+# now fold in data from the other intermediate files
+for filename in other_files:
+    sd = SD(filename)
+    rads_sum += sd.select('MODIS_Radiance_Sum').get()
+    rads_squared_sum += sd.select('MODIS_Radiance_Squared_Sum').get()
+    n += sd.select('MODIS_N').get()
+    flags |= sd.select('MODIS_Flags').get()
+
+# initialize output arrays
+rad_mean = numpy.empty((16, 135, 90), numpy.float64)
+rad_std = numpy.empty((16, 135, 90), numpy.float64)
+bt_mean = numpy.empty((16, 135, 90), numpy.float64)
+std_bt = numpy.empty((16, 135, 90), numpy.float64)
+
+# calculate radiance mean and standard deviation
+rad_mean[:] = rads_sum / n
+rad_std[:] = numpy.sqrt((rads_squared_sum / n) - (rad_mean * rad_mean))
+
+# now convert the radiances to brightness temperature
+bt_mean[:] = modis_bright(rad_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_bt[:] = modis_bright(rad_mean + rad_std) - bt_mean
+
+# for any pixel where N is zero, use fill values for the other
+# variables
+missing_idxs = (n == 0)
+rad_mean[missing_idxs] = -9999.0
+rad_std[missing_idxs] = -9999.0
+bt_mean[missing_idxs] = -9999.0
+std_bt[missing_idxs] = -9999.0
+
+# output to HDF
+
+sd = SD(output_file, SDC.WRITE | SDC.CREATE | SDC.TRUNC)
+
+n_sds = sd.create('MODIS_N', SDC.INT32, n.shape)
+rad_mean_sds = sd.create('MODIS_Radiance', SDC.FLOAT64, rad_mean.shape)
+rad_std_sds = sd.create('MODIS_Radiance_Std', SDC.FLOAT64, rad_std.shape)
+bt_mean_sds = sd.create('MODIS_Brightness_Temp', SDC.FLOAT64, bt_mean.shape)
+std_bt_sds = sd.create('MODIS_Std_Brightness', SDC.FLOAT64, std_bt.shape)
+flags_sds = sd.create('MODIS_Flags', SDC.INT8, flags.shape)
+
+rad_mean_sds.setfillvalue(-9999.0)
+rad_std_sds.setfillvalue(-9999.0)
+bt_mean_sds.setfillvalue(-9999.0)
+std_bt_sds.setfillvalue(-9999.0)
+
+n_sds[:] = n
+rad_mean_sds[:] = rad_mean
+rad_std_sds[:] = rad_std
+bt_mean_sds[:] = bt_mean
+std_bt_sds[:] = std_bt
+flags_sds[:] = flags