"""
pyart.io.sigmet
===============

Reading and writing of Sigmet (raw format) files

.. autosummary::
    :toctree: generated/

    read_sigmet
    ymds_time_to_datetime
    _is_time_ordered_by_reversal
    _is_time_ordered_by_roll
    _is_time_ordered_by_reverse_roll
    _time_order_data_and_metadata_roll
    _time_order_data_and_metadata_reverse
    _time_order_data_and_metadata_full

"""

from __future__ import division
import datetime
import warnings

import numpy as np

from ..config import FileMetadata, get_fillvalue
from ..core.radar import Radar
from .common import make_time_unit_str, _test_arguments, prepare_for_read
from ._sigmetfile import SigmetFile, bin4_to_angle, bin2_to_angle
from . import _sigmet_noaa_hh

SPEED_OF_LIGHT = 299793000.0


def read_sigmet(filename, field_names=None, additional_metadata=None,
                file_field_names=False, exclude_fields=None,
                time_ordered='none', full_xhdr=None, noaa_hh_hdr=None,
                debug=False, ignore_xhdr=False, ignore_sweep_start_ms=None,
                **kwargs):
    """
    Read a Sigmet (IRIS) product file.

    Parameters
    ----------
    filename : str
        Name of Sigmet (IRIS) product file to read or file-like object
        pointing to the beginning of such a file.
    field_names : dict, optional
        Dictionary mapping Sigmet data type names to radar field names. If a
        data type found in the file does not appear in this dictionary or has
        a value of None it will not be placed in the radar.fields dictionary.
        A value of None, the default, will use the mapping defined in the
        metadata configuration file.
    additional_metadata : dict of dicts, optional
        Dictionary of dictionaries to retrieve metadata from during this read.
        This metadata is not used during any successive file reads unless
        explicitly included.  A value of None, the default, will not
        introduct any addition metadata and the file specific or default
        metadata as specified by the metadata configuration file will be used.
    file_field_names : bool, optional
        True to use the Sigmet data type names for the field names. If this
        case the field_names parameter is ignored. The field dictionary will
        likely only have a 'data' key, unless the fields are defined in
        `additional_metadata`.
    exclude_fields : list or None, optional
        List of fields to exclude from the radar object. This is applied
        after the `file_field_names` and `field_names` parameters.
    time_ordered : 'none', 'sequential', 'full', ...,  optional
        Parameter controlling if and how the rays are re-ordered by time.
        The default, 'none' keeps the rays ordered in the same manner as
        they appears in the Sigmet file.  'sequential' will determind and
        apply an operation which maintains a sequential ray order in elevation
        or azimuth yet orders the rays according to time.  If no operation can
        be found to accomplish this a warning is issue and the rays are
        returned in their original order.  'roll', 'reverse', and
        'reverse_and_roll' will apply that operation to the rays in order to
        place them in time order, direct use of these is not recommended.
        'full' will order the rays in strictly time increasing order,
        but the rays will likely become non-sequential, thisoption is not
        recommended unless strict time increasing order is required.
    full_xhdr : bool or None
        Flag to read in all extended headers for possible decoding. None will
        determine if extended headers should be read in automatically by
        examining the extended header type.
    noaa_hh_hdr : bool or None
        Flag indicating if the extended header should be decoded as those
        used by the NOAA Hurricane Hunters aircraft radars.  None will
        determine if the extended header is of this type automatically by
        examining the header. The `full_xhdr` parameter is set to True
        when this parameter is True.
    ignore_xhdr : bool, optional
        True to ignore all data in the extended headers if they exist.
        False, the default, extracts milliseconds precision times and other
        parameter from the extended headers if they exists in the file.
    ignore_sweep_start_ms : bool or None, optional
        True to ignore the millisecond parameter in the start time for each
        sweep, False will uses this parameter when determining the timing of
        each ray. None, the default, will ignore the millisecond sweep start
        timing only when the file does not contain extended headers or when
        the extended header has been explicity ignored using the `ignore_xhdr`
        parameter.  The TRMM RSL library ignores these times so setting this
        parameter to True is required to match the times determined when
        reading Sigmet files with :py:func:`pyart.io.read_rsl`.
        When there are not extended headers ignoring the millisecond sweep
        times provides time data which is always prior to the actual
        collection time with an error from 0 to 2 seconds.
    debug : bool, optional
        Print debug information during read.

    Returns
    -------
    radar : Radar
        Radar object

    """
    # test for non empty kwargs
    _test_arguments(kwargs)

    # create metadata retrieval object
    filemetadata = FileMetadata('sigmet', field_names, additional_metadata,
                                file_field_names, exclude_fields)

    # open the file
    sigmetfile = SigmetFile(prepare_for_read(filename), debug=debug)
    ingest_config = sigmetfile.ingest_header['ingest_configuration']
    task_config = sigmetfile.ingest_header['task_configuration']

    # determine if full extended headers should be read
    if noaa_hh_hdr is True:
        full_xhdr = True
    if full_xhdr is None:
        type_mask = task_config['task_dsp_info']['current_data_type_mask']
        if type_mask['extended_header_type'] == 2:
            full_xhdr = True
        else:
            full_xhdr = False

    # read data
    sigmet_data, sigmet_metadata = sigmetfile.read_data(full_xhdr=full_xhdr)
    first_data_type = sigmetfile.data_type_names[0]
    if first_data_type == 'XHDR':   # don't use XHDR as the first data type
        first_data_type = sigmetfile.data_type_names[1]
    sigmetfile.close()
    nsweeps, nrays, nbins = sigmet_data[first_data_type].shape

    if nsweeps == 0:
        raise IOError('File contains no readable sweep data.')

    # ignore extended header if user requested
    if ignore_xhdr:
        if 'XHDR' in sigmet_data:
            sigmet_data.pop('XHDR')
            sigmet_metadata.pop('XHDR')

    # parse the extended headers for time
    if full_xhdr and 'XHDR' in sigmet_data:
        # extract the ms timing data and store the full header for later
        # analysis, keep both in the sigmet_data dictionary so time ordering
        # and removal of missing rays is performed on these "fields"
        xhdr = sigmet_data.pop('XHDR')
        sigmet_data['XHDR'] = xhdr[:, :, :2].copy().view('i4')
        sigmet_data['XHDR_FULL'] = xhdr
        xhdr_metadata = {}
        for key in sigmet_metadata['XHDR'].keys():
            xhdr_metadata[key] = sigmet_metadata['XHDR'][key].copy()
        sigmet_metadata['XHDR_FULL'] = xhdr_metadata

    # remove missing rays from the data
    good_rays = (sigmet_metadata[first_data_type]['nbins'] != -1)
    rays_missing = (sigmet_metadata[first_data_type]['nbins'] == -1).sum()
    for field_name in sigmet_data.keys():
        sigmet_data[field_name] = sigmet_data[field_name][good_rays]
        field_metadata = sigmet_metadata[field_name]
        for key in field_metadata.keys():
            field_metadata[key] = field_metadata[key][good_rays]
    rays_per_sweep = good_rays.sum(axis=1)

    # time order
    if time_ordered == 'sequential':
        # Determine appropiate time ordering operation for the sweep and
        # determine if that operation will in fact order the times.
        # If it does not issue a warning and perform no time ordering
        if task_config['task_scan_info']['antenna_scan_mode'] == 2:
            # RHI scan
            if _is_time_ordered_by_reversal(
                    sigmet_data, sigmet_metadata, rays_per_sweep):
                time_ordered = 'reverse'
            else:
                warnings.warn('Rays not collected sequentially in time.')
                time_ordered = 'none'
        else:
            # PPI scan
            if _is_time_ordered_by_roll(
                    sigmet_data, sigmet_metadata, rays_per_sweep):
                time_ordered = 'roll'
            elif _is_time_ordered_by_reverse_roll(
                    sigmet_data, sigmet_metadata, rays_per_sweep):
                time_ordered = 'reverse_and_roll'
            else:
                warnings.warn('Rays not collected sequentially in time.')
                time_ordered = 'none'

    if time_ordered == 'full':
        _time_order_data_and_metadata_full(
            sigmet_data, sigmet_metadata, rays_per_sweep)
    if time_ordered == 'reverse':
        _time_order_data_and_metadata_reverse(
            sigmet_data, sigmet_metadata, rays_per_sweep)
    if time_ordered == 'roll':
        _time_order_data_and_metadata_roll(
            sigmet_data, sigmet_metadata, rays_per_sweep)
    if time_ordered == 'reverse_and_roll':
        # reverse followed by roll
        _time_order_data_and_metadata_reverse(
            sigmet_data, sigmet_metadata, rays_per_sweep)
        _time_order_data_and_metadata_roll(
            sigmet_data, sigmet_metadata, rays_per_sweep)

    # sweep_start_ray_index and sweep_end_ray_index
    ray_count = good_rays.sum(axis=1)
    total_rays = ray_count.sum()
    sweep_start_ray_index = filemetadata('sweep_start_ray_index')
    sweep_end_ray_index = filemetadata('sweep_end_ray_index')

    ssri = np.cumsum(np.append([0], ray_count[:-1])).astype('int32')
    sweep_start_ray_index['data'] = ssri
    sweep_end_ray_index['data'] = np.cumsum(ray_count).astype('int32') - 1

    # time
    # The time of each ray collected in a Sigmet file is stored in the ray
    # header as an offset from the start of of the sweep. If the file contains
    # extended ray headers then the time is truncated to millisecond
    # precision. Without extended headers the time is truncated to second
    # precision. The time of the sweep start is recorded in the ingest data
    # header to millisecond percision.
    # The method below sets the volume starting time to the time of the first
    # sweep truncated to the nearest second so that values time['data'] are
    # always positive.  The timing of each ray is then calculated from the
    # offset of the time of the sweep from this volume starting time
    # and the offset from the sweep start to the ray recorded in the ray
    # header.  Additional, the user can select to trucate the sweep start
    # times to second precision using the `ignore_sweep_start_ms` parameter.
    # This may be desired to match libraries such as TRMM RSL, which do not
    # read the millisecond sweep start precision.  By default this
    # trucation is applied when the file contains no extended header.  This
    # ensures that the times listed occur BEFORE the ray was collected with an
    # error from 0 to 2 seconds from the recorded time.  Other methods of
    # combining the mixed percision times would result in times which may have
    # a negative error which is undesired.

    # extract times from the ray headers
    if 'XHDR' in sigmet_data:   # use time in extended headers
        tdata = sigmet_data.pop('XHDR')
        tdata = (tdata.flatten() / 1000.).astype('float64')
        sigmet_extended_header = True
    else:   # use times from the standard ray headers
        tdata = sigmet_metadata[first_data_type]['time'].astype('float64')
        sigmet_extended_header = False

    # determine sweep start times and possibly trucate them to sec. precision
    dts = [ymds_time_to_datetime(d['sweep_start_time'])
           for d in sigmetfile.ingest_data_headers[first_data_type]]
    if ignore_sweep_start_ms is None:
        ignore_sweep_start_ms = not sigmet_extended_header
    if ignore_sweep_start_ms:
        dts = [d.replace(microsecond=0) for d in dts]

    # truncate volume start time to second precision of first sweep
    dt_start = dts[0].replace(microsecond=0)

    # add sweep start times to all time values in each sweep.
    for i, dt in enumerate(dts):
        start = sweep_start_ray_index['data'][i]
        end = sweep_end_ray_index['data'][i]
        td = (dt - dt_start)    # time delta from volume start
        tdata[start:end + 1] += (td.microseconds + (td.seconds + td.days *
                                 24 * 3600) * 10**6) / 10**6

    # stores this data in the time dictionary
    time = filemetadata('time')
    time['data'] = tdata
    time['units'] = make_time_unit_str(dt_start)

    # _range
    _range = filemetadata('range')
    range_info = task_config['task_range_info']
    gate_0 = range_info['first_bin_range'] / 100.       # meters
    gate_nbin = range_info['last_bin_range'] / 100.     # meters
    gate_size = round((gate_nbin - gate_0) / (nbins))
    _range['data'] = gate_0 + gate_size * np.arange(nbins, dtype='float32')
    _range['meters_to_center_of_first_gate'] = np.array([gate_0],
                                                        dtype='float32')
    _range['meters_between_gates'] = np.array([gate_size], dtype='float32')

    # fields
    fields = {}
    for data_type_name, fdata in sigmet_data.items():
        if data_type_name == 'XHDR_FULL':
            continue
        field_name = filemetadata.get_field_name(data_type_name)
        if field_name is None:
            continue
        field_dic = filemetadata(field_name)
        field_dic['data'] = fdata.reshape(-1, nbins)
        field_dic['_FillValue'] = get_fillvalue()
        fields[field_name] = field_dic

    # metadata
    metadata = filemetadata('metadata')
    metadata['original_container'] = 'sigmet'
    metadata['instrument_name'] = ingest_config['site_name'].strip()
    metadata['sigmet_task_name'] = (
        sigmetfile.product_hdr['product_configuration']['task_name'])

    if sigmet_extended_header:
        metadata['sigmet_extended_header'] = 'true'
    else:
        metadata['sigmet_extended_header'] = 'false'
    metadata['time_ordered'] = time_ordered
    metadata['rays_missing'] = rays_missing

    # scan_type
    if task_config['task_scan_info']['antenna_scan_mode'] == 2:
        scan_type = 'rhi'
    else:
        scan_type = 'ppi'

    # latitude
    latitude = filemetadata('latitude')
    lat = bin4_to_angle(ingest_config['latitude_radar'])
    if lat > 180.0:
        lat -= 360.0
    latitude['data'] = np.array([lat], dtype='float64')

    # longitude
    longitude = filemetadata('longitude')
    lon = bin4_to_angle(ingest_config['longitude_radar'])
    if lon > 180.0:
        lon -= 360.0
    longitude['data'] = np.array([lon], dtype='float64')

    # altitude
    altitude = filemetadata('altitude')
    alt = sigmetfile.product_hdr['product_end']['ground_height']
    altitude['data'] = np.array([alt], dtype='float64')

    # sweep_number
    sweep_number = filemetadata('sweep_number')
    sweep_number['data'] = np.arange(nsweeps, dtype='int32')

    # fixed_angle
    fixed_angle = filemetadata('fixed_angle')
    fa = [d['fixed_angle'] for d in
          sigmetfile.ingest_data_headers[first_data_type]]
    fixed_angle['data'] = bin2_to_angle(np.array(fa)).astype('float32')

    # azimuth
    az0 = sigmet_metadata[first_data_type]['azimuth_0']
    az1 = sigmet_metadata[first_data_type]['azimuth_1']
    el0 = sigmet_metadata[first_data_type]['elevation_0']
    el1 = sigmet_metadata[first_data_type]['elevation_1']

    azimuth = filemetadata('azimuth')
    az = (az0 + az1) / 2.
    az[np.where(np.abs(az0 - az1) > 180.0)] += 180.
    az[az > 360.0] -= 360.
    azimuth['data'] = az.astype('float32')

    # elevation
    elevation = filemetadata('elevation')
    elevation['data'] = ((el0 + el1) / 2.).astype('float32')

    # instrument_parameters
    prt = filemetadata('prt')
    prt_mode = filemetadata('prt_mode')
    nyquist_velocity = filemetadata('nyquist_velocity')
    unambiguous_range = filemetadata('unambiguous_range')
    beam_width_h = filemetadata('radar_beam_width_h')
    beam_width_v = filemetadata('radar_beam_width_v')
    pulse_width = filemetadata('pulse_width')

    prt_value = 1. / sigmetfile.product_hdr['product_end']['prf']
    prt['data'] = prt_value * np.ones(total_rays, dtype='float32')

    ur_value = SPEED_OF_LIGHT * prt_value / 2.
    unambiguous_range['data'] = ur_value * np.ones(total_rays, dtype='float32')

    # TODO Multi PRF mode when
    # task_config['task_dsp_info']['multi_prf_flag'] != 0
    prt_mode['data'] = np.array(nsweeps * ['fixed'])

    wavelength_cm = sigmetfile.product_hdr['product_end']['wavelength']
    nv_value = wavelength_cm / (10000.0 * 4.0 * prt_value)
    nyquist_velocity['data'] = nv_value * np.ones(total_rays, dtype='float32')
    beam_width_h['data'] = np.array([bin4_to_angle(
        task_config['task_misc_info']['horizontal_beamwidth'])],
        dtype='float32')
    beam_width_v['data'] = np.array([bin4_to_angle(
        task_config['task_misc_info']['vertical_beamwidth'])],
        dtype='float32')
    pulse_width['data'] = np.array(
        [task_config['task_dsp_info']['pulse_width'] * 1e-8] *
        len(time['data']), dtype='float32')

    instrument_parameters = {'unambiguous_range': unambiguous_range,
                             'prt_mode': prt_mode, 'prt': prt,
                             'nyquist_velocity': nyquist_velocity,
                             'radar_beam_width_h': beam_width_h,
                             'radar_beam_width_v': beam_width_v,
                             'pulse_width': pulse_width}

    # decode extended headers
    extended_header_params = {}
    if noaa_hh_hdr is None:
        type_mask = task_config['task_dsp_info']['current_data_type_mask']
        htype = type_mask['extended_header_type']
        if full_xhdr and htype == 2 and sigmet_data['XHDR_FULL'][0, 3] == 136:
            noaa_hh_hdr = True
        else:
            noaa_hh_hdr = False

    if noaa_hh_hdr:
        t = _sigmet_noaa_hh._decode_noaa_hh_hdr(
            sigmet_data['XHDR_FULL'], filemetadata, azimuth, elevation)
        (latitude, longitude, altitude, extended_header_params) = t
        metadata['platform_type'] = 'aircraft'
        # scan_type determined from the antenna_scan_mode parameters
        noaa_hh_scan_modes = {4: 'ppi', 7: 'rhi'}
        scan_mode = task_config['task_scan_info']['antenna_scan_mode']
        if scan_mode not in noaa_hh_scan_modes:
            warnings.warn("Unknown antenna_scan_mode, defaulting to 'rhi'.")
            scan_type = 'rhi'
        else:
            scan_type = noaa_hh_scan_modes[scan_mode]

    # sweep_mode
    sweep_mode = filemetadata('sweep_mode')
    if scan_type == 'ppi':
        sweep_mode['data'] = np.array(nsweeps * ['azimuth_surveillance'])
    else:
        sweep_mode['data'] = np.array(nsweeps * ['rhi'])

    return Radar(
        time, _range, fields, metadata, scan_type,
        latitude, longitude, altitude,
        sweep_number, sweep_mode, fixed_angle, sweep_start_ray_index,
        sweep_end_ray_index,
        azimuth, elevation,
        instrument_parameters=instrument_parameters,
        **extended_header_params)


def _is_time_ordered_by_reversal(data, metadata, rays_per_sweep):
    """
    Returns if volume can be time ordered by reversing some or all sweeps.
    True if the volume can be time ordered, False if not.
    """
    if 'XHDR' in data:
        ref_time = data['XHDR'].copy()
        ref_time.shape = ref_time.shape[:-1]
    else:
        ref_time = metadata[metadata.keys()[0]]['time'].astype('int32')
    start = 0
    for nrays in rays_per_sweep:
        if nrays == 0 or nrays == 1:
            continue    # Do not attempt to order sweeps with no rays
        s = slice(start, start + nrays)     # slice which selects sweep
        start += nrays
        sweep_time_diff = np.diff(ref_time[s])
        if np.all(sweep_time_diff >= 0) or np.all(sweep_time_diff <= 0):
            continue
        else:
            return False
    return True


def _is_time_ordered_by_roll(data, metadata, rays_per_sweep):
    """
    Returns if volume can be time ordered by rolling some or all sweeps.
    True if the volume can be time ordered, False if not.
    """
    if 'XHDR' in data:
        ref_time = data['XHDR'].copy()
        ref_time.shape = ref_time.shape[:-1]
    else:
        ref_time = metadata[metadata.keys()[0]]['time'].astype('int32')
    start = 0
    for nrays in rays_per_sweep:
        if nrays == 0 or nrays == 1:
            continue    # Do not attempt to order sweeps with no rays
        s = slice(start, start + nrays)     # slice which selects sweep
        first = ref_time[start]
        last = ref_time[start + nrays - 1]
        start += nrays
        sweep_time_diff = np.diff(ref_time[s])
        count = np.count_nonzero(sweep_time_diff < 0)
        # compare the first and last times for continuity
        if (first - last) < 0:
            count += 1
        if count != 0 and count != 1:
            return False
    return True


def _is_time_ordered_by_reverse_roll(data, metadata, rays_per_sweep):
    """
    Returns if volume can be time ordered by reversing and rolling some or all
    sweeps.  True if the volume can be time ordered, False if not.
    """
    if 'XHDR' in data:
        ref_time = data['XHDR'].copy()
        ref_time.shape = ref_time.shape[:-1]
    else:
        ref_time = metadata[metadata.keys()[0]]['time'].astype('int32')
    start = 0
    for nrays in rays_per_sweep:
        if nrays == 0 or nrays == 1:
            continue    # Do not attempt to order sweeps with no rays
        s = slice(start, start + nrays)     # slice which selects sweep
        first = ref_time[start]
        last = ref_time[start + nrays - 1]
        start += nrays
        sweep_time_diff = np.diff(ref_time[s])
        if sweep_time_diff.min() < 0:   # optional reverse
            sweep_time_diff = np.diff(ref_time[s][::-1])
            first, last = last, first
        count = np.count_nonzero(sweep_time_diff < 0)
        # compare the first and last times for continuity
        if (first - last) < 0:
            count += 1
        if count != 0 and count != 1:
            return False
    return True


def _time_order_data_and_metadata_roll(data, metadata, rays_per_sweep):
    """
    Put Sigmet data and metadata in time increasing order using a roll
    operation.
    """
    # Sigmet data is stored by sweep in azimuth or elevation increasing order.
    # Time ordering PPI scans can typically be achieved by rolling the
    # ray collected first to the beginning of the sweep, which is performed
    # here.  Perfect time ordering is achieved if the rays within the sweep
    # were collected sequentially in a clockwise manner from 0 to 360 degrees
    # regardless of the first azimuth collected.
    if 'XHDR' in data:
        ref_time = data['XHDR'].copy()
        ref_time.shape = ref_time.shape[:-1]
    else:
        ref_time = metadata[metadata.keys()[0]]['time'].astype('int32')

    start = 0
    for nrays in rays_per_sweep:
        if nrays == 0 or nrays == 1:
            continue    # Do not attempt to order sweeps with no rays

        s = slice(start, start + nrays)     # slice which selects sweep
        start += nrays
        # determine the number of place by which elements should be shifted.
        sweep_time = ref_time[s]
        sweep_time_diff = np.diff(sweep_time)
        if sweep_time_diff.min() >= 0:
            continue    # already time ordered
        shift = -(sweep_time_diff.argmin() + 1)

        # roll the data and metadata for each field
        for field in data.keys():
            data[field][s] = np.roll(data[field][s], shift, axis=0)
            field_metadata = metadata[field]
            for key in field_metadata.keys():
                field_metadata[key][s] = np.roll(field_metadata[key][s], shift)
    return


def _time_order_data_and_metadata_reverse(data, metadata, rays_per_sweep):
    """
    Put Sigmet data and metadata in time increasing order by reverse sweep in
    time reversed order.
    """
    # Sigmet data is stored by sweep in azimuth or elevation increasing order.
    # Time ordering RHI scans can typically be achieved by reversing the
    # ray order of sweep collected in time from 180 to 0 degrees.
    # Perfect time ordering is achieved if the rays within all sweeps
    # were collected sequentially from 0 to 180 degree or 180 to 0 degrees.
    if 'XHDR' in data:
        ref_time = data['XHDR'].copy()
        ref_time.shape = ref_time.shape[:-1]
    else:
        ref_time = metadata[metadata.keys()[0]]['time'].astype('int32')

    start = 0
    for nrays in rays_per_sweep:
        if nrays == 0 or nrays == 1:
            continue    # Do not attempt to order sweeps with to few rays

        s = slice(start, start + nrays)     # slice which selects sweep
        start += nrays
        # determine the number of place by which elements should be shifted.
        sweep_time = ref_time[s]
        sweep_time_diff = np.diff(sweep_time)
        if sweep_time_diff.min() >= 0:
            continue    # already time ordered, no reversal needed
        # reverse the data and metadata for each field
        for field in data.keys():
            data[field][s] = data[field][s][::-1]
            field_metadata = metadata[field]
            for key in field_metadata.keys():
                field_metadata[key][s] = field_metadata[key][s][::-1]
    return


def _time_order_data_and_metadata_full(data, metadata, rays_per_sweep):
    """
    Put Sigmet data and metadata in time increasing order by sorting the
    times.
    """
    # Sigmet data is stored by sweep in azimuth or elevation increasing order.
    # When rays within the sweeps are collected non-sequentially or in a
    # complex manner, perfect time ordering can only be achived by sorting
    # the times themselves and reordering the rays according to this sort.
    # This ordering method should only be used as a last resort when perfect
    # time ordering is required in the output and other ordering operations
    # (roll, reverse, reverse-roll) will not order the rays correctly.
    if 'XHDR' in data:
        ref_time = data['XHDR'].copy()
        ref_time.shape = ref_time.shape[:-1]
    else:
        ref_time = metadata[metadata.keys()[0]]['time'].astype('int32')

    start = 0
    for nrays in rays_per_sweep:
        if nrays == 0 or nrays == 1:
            continue    # Do not attempt to order sweeps with no rays

        s = slice(start, start + nrays)     # slice which selects sweep
        start += nrays
        # determine the indices which sort the sweep time using a stable
        # sorting algorithm to prevent excessive azimuth scrambling.
        sweep_time = ref_time[s]
        sweep_time_diff = np.diff(ref_time[s])
        if sweep_time_diff.min() >= 0:
            continue    # already time ordered
        sort_idx = np.argsort(sweep_time, kind='mergesort')

        # sort the data and metadata for each field
        for field in data.keys():
            data[field][s] = data[field][s][sort_idx]
            field_metadata = metadata[field]
            for key in field_metadata.keys():
                field_metadata[key][s] = field_metadata[key][s][sort_idx]
    return


def ymds_time_to_datetime(ymds):
    """ Return a datetime object from a Sigmet ymds_time dictionary. """
    dt = datetime.datetime(ymds['year'], ymds['month'], ymds['day'])
    # lowest 10 bits of millisecond parameter specifies milliseconds
    microsec = 1e3 * (ymds['milliseconds'] & 0b1111111111)
    delta = datetime.timedelta(seconds=ymds['seconds'], microseconds=microsec)
    return dt + delta
