Add mirror / detector separation to iasi_modis

intercal/iasi_modis.py

@@ -161,6 +161,9 @@ def sge_main(case, input_dir):
 
     result_names = ['Latitude', 'Longitude', 'Time',
                     'N', 'Mean_T_b', 'T_b_Std']
+    if sensor_cls is Modis:
+        result_names += ['Mirror_Detector_N', 'Mirror_Detector_Mean_T_b',
+                         'Mirror_Detector_T_b_Std']
     results = {k: [] for k in result_names}
     def result(name, value):
         results[name].append(value)
@@ -188,6 +191,23 @@ def sge_main(case, input_dir):
         result('Mean_T_b', mu)
         result('T_b_Std', sigma)
 
+        if isinstance(sensor, Modis):
+            md_n = np.empty([2, 10], np.int)
+            md_mu = np.empty([2, 10], np.float)
+            md_sigma = np.empty([2, 10], np.float)
+            for mirror in range(2):
+                for detector in range(10):
+                    mask = sensor.detector_mask(detector)
+                    mask[modis.mirror_side != mirror] = False
+                    n, mu, sigma = process_sno(sensor, sno.lat, sno.lon,
+                                               bands, d_max, mask)
+                    md_n[mirror,detector] = n
+                    md_mu[mirror,detector] = mu
+                    md_sigma[mirror,detector] = sigma
+            result('Mirror_Detector_N', md_n)
+            result('Mirror_Detector_Mean_T_b', md_mu)
+            result('Mirror_Detector_T_b_Std', md_sigma)
+
         print 'DONE', sno_str(sno)
 
     if not results['Latitude']:
@@ -199,12 +219,19 @@ def sge_main(case, input_dir):
             out_file.create_dataset(name, data=np.array(results[name]))
     print 'RESULTS IN', out_filename
 
-def process_sno(sensor, lat, lon, bands, d_max):
+def process_sno(sensor, lat, lon, bands, d_max, pre_mask=None):
+
+    sensor_lat = sensor.latitude
+    sensor_lon = sensor.longitude
+    if pre_mask:
+        sensor_lat = sensor_lat[pre_mask]
+        sensor_lon = sensor_lon[pre_mask]
 
-    d = gcdist(sensor.latitude, sensor.longitude, lat, lon)
+    d = gcdist(sensor_lat, sensor_lon, lat, lon)
     mask = d <= d_max
 
     if isinstance(sensor, Iasi):
+        assert not pre_mask
         R = sensor.radiance()[mask]
 
     n, mu, sigma = [], [], []
@@ -218,7 +245,10 @@ def process_sno(sensor, lat, lon, bands, d_max):
             R_band = srf.apply(sensor.wavenumber, R)
 
         else:
-            R_band = sensor.radiance(band)[mask]
+            R = sensor.radiance(band)
+            if pre_mask:
+                R = R[pre_mask]
+            R_band = R[mask]
 
         # filter out NaN fill values
         R_band = R_band[np.logical_not(np.isnan(R_band))]

intercal/modis.py

@@ -5,7 +5,7 @@ from datetime import datetime, timedelta
 
 import numpy as np
 from pkg_resources import resource_filename
-from pyhdf import SD
+from pyhdf import HDF, SD, VS
 
 from .sensor import Sensor
 from .srf import WavelengthSrf
@@ -82,6 +82,30 @@ class Modis(Sensor):
         rad[mask] = np.nan
         sd.end()
 
+    def detector_mask(self, detector):
+
+        assert detector in range(10)
+        mask = np.zeros([self.shape[0]], np.bool)
+        mask[detector::10] = True
+        return mask
+
+    @property
+    def mirror_side(self):
+
+        if not hasattr(self, 'mirror'):
+
+            # FIXME: how to properly close the HDF handle?
+            mirror_sides = []
+            for f in self.files:
+                vs = HDF.HDF(f['MOD021KM']).vstart()
+                vd = vs.attach('Level 1B Swath Metadata')
+                vd.setfields('Mirror Side')
+                mirror_sides.extend(a[0] for a in vd[:])
+
+            self.mirror = np.repeat(np.array(mirror_sides), 10)
+
+        return self.mirror
+
 def modis_srf(band):
 
     """Return the Terra MODIS SRF for the given emissive band"""