Add h5cat function

intercal/iasi_modis.py

@@ -14,6 +14,45 @@ from .sensor import product_subdir
 from .snotimes import SnoTimes
 from .util import gcdist
 
+bands = [20, 21, 22, 23, 24, 25, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]
+d_max = 50.0
+
+def h5cat(dir, pattern, out_filename):
+
+    """Concatenate partial HDF5 result files into one"""
+
+    in_filenames = glob(os.path.join(dir, pattern))
+
+    # find the total size of the first dimension across all files
+    # (we assume the first dimension of each dataset is the same,
+    # which it should be)
+    size = 0
+    for in_filename in in_filenames:
+        with h5py.File(in_filename, 'r') as in_file:
+            size += len(in_file['Latitude'])
+
+    # output file goes in same directory with inputs
+    with h5py.File(os.path.join(dir, out_filename), 'w') as out_file:
+
+        # use first input file to determine the datasets to copy (this
+        # assumes there are no groups)
+        out_datasets = {}
+        with h5py.File(in_filenames[0], 'r') as in_file:
+            idx_fin = len(in_file['Latitude'])
+            for name, in_dataset in in_file.items():
+                shape = [size] + list(in_dataset.shape[1:])
+                out_datasets[name] = out_file.create_dataset(
+                                         name, shape, in_dataset.dtype)
+                out_datasets[name][:idx_fin] = in_dataset
+
+        # now copy over the rest of the files as well
+        for in_filename in in_filenames[1:]:
+            with h5py.File(in_filename, 'r') as in_file:
+                idx_ini = idx_fin
+                idx_fin = idx_ini + len(in_file['Latitude'])
+                for name, out_dataset in out_datasets.items():
+                    out_dataset[idx_ini:idx_fin] = in_file[name]
+
 def sge_main(case, input_dir):
 
     """Driver function to be called from an SGE job"""
@@ -43,9 +82,6 @@ def sge_main(case, input_dir):
     assert snotimes.platform_1 == 'METOPA'
     assert snotimes.platform_2 == 'TERRA'
 
-    bands = [20, 21, 22, 23, 24, 25, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]
-    d_max = 50.0
-
     result_names = ['Latitude', 'Longitude', 'Time',
                     'N', 'Mean_T_b', 'T_b_Std']
     results = {k: [] for k in result_names}