Add crisiasiairs package

src/intercal/crisiasiairs/airs.py

+
+import numpy as np
+import pyhdf.HDF
+import pyhdf.SD
+import pyhdf.VS
+
+class AirsGranule(object):
+
+    num_fors_per_scan = 90
+    num_fovs_per_for = 1
+
+    has_imaginary_radiance = False
+
+    num_files = 1
+
+    def __init__(self, filename):
+        self.sd = pyhdf.SD.SD(filename)
+        self.wavenumber = self.read_wavenumber(filename)
+
+    def read_wavenumber(self, filename):
+        num_channels = 2378
+        nested_list = pyhdf.HDF.HDF(filename).vstart().attach('nominal_freq').read(num_channels)
+        return np.array(nested_list).flatten()
+
+    @property
+    def num_scans(self):
+        return self.sd.datasets()['radiances'][1][0]
+
+    def read_latitude(self, latitude):
+        latitude[:,:,0] = self.sd.select('Latitude')[:]
+
+    def read_longitude(self, longitude):
+        longitude[:,:,0] = self.sd.select('Longitude')[:]
+
+    def read_radiance(self, radiance, imaginary):
+        if imaginary:
+            raise ValueError('no imaginary radiances for AIRS')
+        radiance[:,:,0,:] = self.sd.select('radiances')[:]
+

src/intercal/crisiasiairs/cris.py

+
+import h5py
+import numpy as np
+
+class CrisGranule(object):
+
+    num_fors_per_scan = 30
+    num_fovs_per_for = 9
+
+    wavenumber = np.concatenate([np.linspace(648.75, 1096.25, 717),
+                                 np.linspace(1207.5, 1752.5, 437),
+                                 np.linspace(2150, 2555, 163)])
+
+    has_imaginary_radiance = True
+
+    num_files = 2
+
+    def __init__(self, gcrso_file, scris_file):
+
+        self.gcrso_file = gcrso_file
+        self.scris_file = scris_file
+
+    @property
+    def num_scans(self):
+
+        with h5py.File(self.gcrso_file, 'r') as gcrso:
+            return len(gcrso['/All_Data/CrIS-SDR-GEO_All/Latitude'])
+
+    def read_latitude(self, latitude):
+
+        with h5py.File(self.gcrso_file, 'r') as gcrso:
+            latitude[:,:,:] = gcrso['/All_Data/CrIS-SDR-GEO_All/Latitude'][:,:,:]
+
+    def read_longitude(self, longitude):
+
+        with h5py.File(self.gcrso_file, 'r') as gcrso:
+            longitude[:,:,:] = gcrso['/All_Data/CrIS-SDR-GEO_All/Longitude'][:,:,:]
+
+    def read_radiance(self, radiance, imaginary):
+
+        with h5py.File(self.scris_file, 'r') as scris:
+            group = scris['/All_Data/CrIS-SDR_All']
+            wavenumber_idx_ini = 0
+            for band in ['LW', 'MW', 'SW']:
+                dataset = 'ES_' + {False: 'Real', True: 'Imaginary'}[imaginary] + band
+                wavenumber_idx_fin = wavenumber_idx_ini + group[dataset].shape[-1]
+                radiance[:,:,:,wavenumber_idx_ini:wavenumber_idx_fin] = group[dataset]
+                wavenumber_idx_ini = wavenumber_idx_fin
+        

src/intercal/crisiasiairs/iasi.py

+
+import functools
+
+import coda
+import numpy as np
+
+class IasiGranule(object):
+
+    wavenumber = np.linspace(645.0, 2760.0, 8461)
+    num_fors_per_scan = 30
+    num_fovs_per_for = 4
+
+    has_imaginary_radiance = False
+
+    num_files = 1
+
+    def __init__(self, filename):
+
+        self.filename = filename
+
+    @property
+    def num_scans(self):
+
+        with CodaHandle(self.filename) as handle:
+            num_scans, = handle.get_size('MDR')
+        return num_scans
+
+    def read_latitude(self, latitude):
+
+        self._read_geolocation('latitude', latitude)
+
+    def read_longitude(self, longitude):
+
+        self._read_geolocation('longitude', longitude)
+
+    def read_radiance(self, radiance, imaginary):
+
+        if imaginary:
+            raise ValueError('Imaginary radiance not available for IASI')
+        with CodaHandle(self.filename) as handle:
+            scaled_radiances = (
+                np.array([a for a in handle.fetch('MDR', -1, 'MDR', 'GS1cSpect')]))
+            scale_factors = handle.fetch('GIADR_ScaleFactors')
+            quality_flags = (
+                np.array([a for a in handle.fetch('MDR', -1, 'MDR', 'GQisFlagQual')],
+                         np.bool))
+        self._unscale_radiances(scaled_radiances, scale_factors, radiance)
+        self._set_low_quality_spectra_to_nan(quality_flags, radiance)
+
+    def _read_geolocation(self, variable, output_array):
+
+        latlon_idx = dict(latitude=1, longitude=0)[variable]
+        with CodaHandle(self.filename) as handle:
+            output_array[:,:,:] = [
+                a for a in handle.fetch('MDR', -1, 'MDR', 'GGeoSondLoc',
+                                        [-1, -1, latlon_idx])]
+
+    def _unscale_radiances(self, scaled_radiances, scale_factors, radiances):
+
+        num_scale_factors = scale_factors.IDefScaleSondNbScale
+        for scale_factor_idx in range(num_scale_factors):
+            scale_factor = scale_factors.IDefScaleSondScaleFactor[scale_factor_idx]
+            wavenumber_idx_ini = (scale_factors.IDefScaleSondNsfirst[scale_factor_idx] -
+                                      self._first_channel_number)
+            wavenumber_idx_fin = (scale_factors.IDefScaleSondNslast[scale_factor_idx] -
+                                      self._first_channel_number + 1)
+            radiances[:,:,:,wavenumber_idx_ini:wavenumber_idx_fin] = (
+                scaled_radiances[:,:,:,wavenumber_idx_ini:wavenumber_idx_fin] *
+                    (10.0 ** -scale_factor * self._radiance_conversion_factor))
+
+    def _set_low_quality_spectra_to_nan(self, quality_flags, radiances):
+
+        radiances[quality_flags.any(axis=-1)] = np.nan
+
+    @property
+    def _first_channel_number(self):
+
+        with CodaHandle(self.filename) as handle:
+            return handle.fetch('MDR', 0, 'MDR', 'IDefNsfirst1b')
+
+    # to convert from W / m^2 / sr^1 / m^-1 to mW / m^2 / sr / cm^-1
+    _radiance_conversion_factor = 1e5
+
+class CodaHandle(object):
+
+    def __init__(self, filename):
+
+        self.handle = coda.open(filename)
+
+    def __enter__(self):
+
+        return self
+
+    def __exit__(self, exc_type, exc_value, exc_trace):
+
+        coda.close(self.handle)
+
+    def __getattr__(self, attr):
+
+        return functools.partial(getattr(coda, attr), self.handle)
+

src/intercal/crisiasiairs/main.py

+
+import argparse
+from datetime import datetime
+import os
+
+import h5py
+import numpy as np
+
+from .airs import AirsGranule
+from .cris import CrisGranule
+from .iasi import IasiGranule
+from .super_granule import SuperGranule
+from .util import gcdist
+
+def _parse_args():
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('platform_1')
+    parser.add_argument('time_1', type=parse_datetime)
+    parser.add_argument('platform_2')
+    parser.add_argument('time_2', type=parse_datetime)
+    parser.add_argument('latitude', type=float)
+    parser.add_argument('longitude', type=float)
+    parser.add_argument('max_distance', type=float)
+    parser.add_argument('sensor', choices=['airs', 'cris', 'iasi'])
+    parser.add_argument('files', nargs='+', metavar='file')
+    return parser.parse_args()
+
+def parse_datetime(s):
+
+    return datetime.strptime(s, '%Y-%m-%dT%H:%M:%S.%f')
+
+def _create_granule(sensor, files):
+
+    granule_cls = {'airs': AirsGranule, 'cris': CrisGranule, 'iasi': IasiGranule}[sensor]
+    return SuperGranule(granule_cls, files)
+
+def _filter_by_proximity(granule, lat, lon, max_distance):
+
+    distance = gcdist(granule.latitude, granule.longitude, lat, lon)
+    return (distance <= max_distance).nonzero()
+
+def _process_spectra(granule, radiance, scan_idx, for_idx, fov_idx,
+                     output_file, dataset_prefix=None):
+
+    radiance = radiance.astype(np.float64)
+    R = _grab_spectra(granule, radiance, scan_idx, for_idx, fov_idx)
+    n, sum, sum_of_squares = _aggregate_by_fov(granule, fov_idx, R)
+    new_dataset_prefix = (
+        'Radiance' if dataset_prefix is None else dataset_prefix + '_Radiance')
+    _output(n, sum, sum_of_squares, output_file, 'FOV_' + new_dataset_prefix)
+    _output(n.sum(axis=0), sum.sum(axis=0), sum_of_squares.sum(axis=0),
+            output_file, new_dataset_prefix)
+
+def _grab_spectra(granule, radiance, scan_idx, for_idx, fov_idx):
+
+    flat_fov_idx = ((scan_idx * granule.num_fors_per_scan +
+                        for_idx) * granule.num_fovs_per_for +
+                        fov_idx)
+    wn_idx = np.arange(granule.num_wavenumbers)
+    flat_sample_idx = (flat_fov_idx[:,None] * granule.num_wavenumbers + wn_idx).flatten()
+    return radiance.take(flat_sample_idx)
+
+def _aggregate_by_fov(granule, fov_idx, radiance):
+
+    wn_idx = np.arange(granule.num_wavenumbers)
+    flat_output_idx = (fov_idx[:,None] * granule.num_wavenumbers + wn_idx).flatten()
+    num_output_idxs = granule.num_fovs_per_for * granule.num_wavenumbers
+    n = np.bincount(fov_idx, minlength=granule.num_fovs_per_for)
+    r = (np.bincount(flat_output_idx, weights=radiance, minlength=num_output_idxs)
+             .reshape([granule.num_fovs_per_for, granule.num_wavenumbers]))
+    r2 = (np.bincount(flat_output_idx, weights=(radiance * radiance),
+                      minlength=num_output_idxs)
+              .reshape([granule.num_fovs_per_for, granule.num_wavenumbers]))
+    return n, r, r2
+
+def _output(n, sum, sum_of_squares, output_file, dataset_prefix):
+
+    mean = sum / n[...,None]
+    std = np.sqrt(sum_of_squares / n[...,None] - mean * mean)
+    output_file.create_dataset(dataset_prefix + '_Count', data=n)
+    output_file.create_dataset(dataset_prefix + '_Mean', data=mean)
+    output_file.create_dataset(dataset_prefix + '_Std', data=std)
+
+def unix_time(t):
+
+    return (t - datetime(1970, 1, 1)).total_seconds()
+
+def output_file_name(args):
+
+    platform_1, platform_2 = sorted([args.platform_1, args.platform_2])
+    earlier_time = min(args.time_1, args.time_2)
+    return 'sno_{}_{}_{}_{}.h5'.format(platform_tag(platform_1).lower(),
+                                       platform_tag(platform_2).lower(),
+                                       args.sensor.lower(),
+                                       earlier_time.strftime('%Y-%m-%dT%H-%M'))
+
+def platform_tag(platform):
+    tags = {'AQUA': 'Aqua', 'METOP-A': 'Metop-A', 'METOP-B': 'Metop-B', 'SUOMI NPP': 'NPP'}
+    return tags[platform]
+
+def main():
+
+    args = _parse_args()
+    granule = _create_granule(args.sensor, args.files)
+    scan_idx, for_idx, fov_idx = (
+        _filter_by_proximity(granule, args.latitude, args.longitude, args.max_distance))
+    with h5py.File(output_file_name(args), 'w') as output_file:
+        _process_spectra(granule, granule.radiance, scan_idx, for_idx, fov_idx,
+                         output_file)
+        if granule.has_imaginary_radiance:
+            _process_spectra(
+                granule, granule.imaginary_radiance, scan_idx, for_idx, fov_idx,
+                output_file, 'Imaginary')
+        output_file.create_dataset('{}_Time'.format(platform_tag(args.platform_1)),
+                                   data=unix_time(args.time_1))
+        output_file.create_dataset('{}_Time'.format(platform_tag(args.platform_2)),
+                                   data=unix_time(args.time_2))
+        output_file.create_dataset('SNO_Latitude', data=args.latitude)
+        output_file.create_dataset('SNO_Longitude', data=args.longitude)
+        output_file.create_dataset('Input_Files', data=[os.path.basename(f) for f in args.files])
+
+if __name__ == '__main__':
+    main()
+

src/intercal/crisiasiairs/super_granule.py

+
+import numpy as np
+
+class SuperGranule(object):
+
+    def __init__(self, granule_cls, files):
+
+        self.granules = [granule_cls(*args)
+                             for args
+                             in self._grouped(files, granule_cls.num_files)]
+
+    @property
+    def num_scans(self):
+
+        return sum(g.num_scans for g in self.granules)
+
+    @property
+    def num_fors_per_scan(self):
+
+        return self.granules[0].num_fors_per_scan
+
+    @property
+    def num_fovs_per_for(self):
+
+        return self.granules[0].num_fovs_per_for
+
+    @property
+    def num_wavenumbers(self):
+
+        return len(self.wavenumber)
+
+    @property
+    def wavenumber(self):
+
+        return self.granules[0].wavenumber
+
+    @property
+    def latitude(self):
+
+        if not hasattr(self, '_latitude'):
+            self._latitude = self._get_geolocation('latitude')
+        return self._latitude
+
+    @property
+    def longitude(self):
+
+        if not hasattr(self, '_longitude'):
+            self._longitude = self._get_geolocation('longitude')
+        return self._longitude
+
+    @property
+    def radiance(self):
+
+        if not hasattr(self, '_radiance'):
+            self._radiance = self._get_radiance()
+        return self._radiance
+
+    @property
+    def has_imaginary_radiance(self):
+
+        return self.granules[0].has_imaginary_radiance
+
+    @property
+    def imaginary_radiance(self):
+
+        if not hasattr(self, '_imaginary_radiance'):
+            self._imaginary_radiance = self._get_radiance(imaginary=True)
+        return self._imaginary_radiance
+
+    def _get_geolocation(self, variable):
+
+        output_array = (
+            np.empty([self.num_scans, self.num_fors_per_scan, self.num_fovs_per_for],
+                     np.float32))
+        for granule, scan_idx_ini, scan_idx_fin in self._granules_with_indexes:
+            read_method = getattr(granule, 'read_' + variable)
+            read_method(output_array[scan_idx_ini:scan_idx_fin,:,:])
+        return output_array
+
+    def _get_radiance(self, imaginary=False):
+
+        output_array = np.empty([self.num_scans, self.num_fors_per_scan,
+                                 self.num_fovs_per_for, len(self.wavenumber)],
+                                np.float32)
+        for granule, scan_idx_ini, scan_idx_fin in self._granules_with_indexes:
+            granule.read_radiance(output_array[scan_idx_ini:scan_idx_fin,:,:],
+                                  imaginary)
+        return output_array
+
+    @property
+    def _granules_with_indexes(self):
+
+        scan_idx = 0
+        for granule in self.granules:
+            yield granule, scan_idx, scan_idx + granule.num_scans
+            scan_idx += granule.num_scans
+
+    def _grouped(self, iterable, n):
+
+        return zip(*([iter(iterable)] * n))
+

src/intercal/crisiasiairs/test/test_iasi.py

+
+from unittest import TestCase
+
+from intercal.iasi import Iasi, CodaFile
+
+class IasiTest(TestCase):
+
+    def test_latitude(self):
+
+        lat = Iasi(EXAMPLE_IASI_FILE).latitude
+        self.assertEqual(lat.shape, (22, 30, 4))
+        self.assertAlmostEqual(lat[0,0,0], 150.069, places=3)
+
+    def test_longitude(self):
+
+        lon = Iasi(EXAMPLE_IASI_FILE).longitude
+        self.assertEqual(lon.shape, (22, 30, 4))
+        self.assertAlmostEqual(lon[0,0,0], -50.907, places=3)
+
+    def test_radiance(self):
+
+        rad = Iasi(EXAMPLE_IASI_FILE).radiance
+
+class CodaFileTest(TestCase):
+
+    def test_getattr(self):
+
+        with CodaFile(EXAMPLE_IASI_FILE) as coda_file:
+            self.assertEqual(coda_file.get_product_filename(), EXAMPLE_IASI_FILE)
+
+EXAMPLE_IASI_FILE = ('/data/gregq/iasi/IASI_xxx_1C_M02_20121106234457Z_' +
+                         '20121106234753Z_N_O_20121107004428Z__20121107004545')
+

src/intercal/crisiasiairs/util.py

+
+"""Utility functions for inter-calibration"""
+
+import numpy as np
+
+# mean earth radius; km
+r_earth_mean = 6371.0
+
+# speed of light; m * s^-1
+c = 2.99792458e8
+
+# Planck constant; J * s
+h = 6.62606957e-34
+
+# Boltzmann constant; J * K^-1
+k_B = 1.3806488e-23
+
+# Planck function coefficients
+c1 = 2.0 * h * c**2
+c2 = h * c / k_B
+
+def gcdist(lat_a, lon_a, lat_b, lon_b):
+
+    """Great circle distance in km between lat lon pair"""
+
+    lat_a, lon_a, lat_b, lon_b = (
+        np.radians(x) for x in [lat_a, lon_a, lat_b, lon_b])
+    tmp = (np.sin(lat_a) * np.sin(lat_b) +
+               np.cos(lat_a) * np.cos(lat_b) * np.cos(lon_a - lon_b))
+    tmp[tmp > 1.0] = 1.0
+    return np.arccos(tmp) * r_earth_mean
+
+def wl2wn(lambda_):
+
+    """Convert from wavelength to wavenumber.
+
+    lambda_ given in microns. Result in inverse centimeters.
+    """
+
+    return 1e4 / lambda_
+
+def wn2wl(nu):
+
+    """Convert from wavenumber to wavelength.
+
+    nu given in inverse centimeters. Result in microns.
+    """
+
+    return 1e4 / nu
+
+def planck_wl(lambda_, T):
+
+    """Planck function for lambda_ in microns, T in Kelvin
+
+    Output in W per m^2 per sr per micron.
+    """
+
+    lambda_ = lambda_ * 1e-6
+    return c1 / lambda_**5 / (np.exp(c2 / (lambda_ * T)) - 1.0) * 1e-6
+
+def planck_wn(nu, T):
+
+    """Planck function for nu in inverse centimeters, T in Kelvin
+
+    Output in mW per m^2 per sr per inverse centimeter.
+    """
+
+    nu = nu * 1e2
+    return c1 * nu**3 / (np.exp(c2 * nu / T) - 1.0) * (1e2 * 1e3)
+
+def inv_planck_wl(lambda_, B):
+
+    """Inverse planck routine for wavelength in microns
+
+    B is radiance in W per m^2 per sr per micron.
+    """
+
+    lambda_ = lambda_ * 1e-6
+    B = B * 1e6
+    return c2 / lambda_ / np.log(c1 / (lambda_**5 * B) + 1.0)
+
+def inv_planck_wn(nu, B):
+
+    """Inverse planck routine for wavenumber in inverse centimeters
+
+    B is radiance in mW per m^2 per sr per inverse centimeter.
+    """
+
+    nu = nu * 1e2
+    B = B * (1e-2 * 1e-3)
+    return c2 * nu / np.log(c1 * nu**3 / B + 1.0)
+
+def rad_wl2wn(lambda_, R):
+
+    """Convert radiances from per-wavelength to per-wavenumber units
+
+    lambda_ is in microns and R is in W per m^2 per sr per micron.
+    Output is in mW per m^2 per sr per inverse centimeter.
+    """
+
+    return R * lambda_**2 * (1e3 * 1e-4)
+
+def rad_wn2wl(nu, R):
+
+    """Convert radiances from per-wavenumber to per-wavelength units
+
+    nu is in inverse centimeters and R is in mW per m^2 per sr per
+    inverse centimeter. Output is in W per m^2 per sr per micron.
+    """
+
+    return R * nu**2 * (1e-3 * 1e-4)
+