class MatchupError(Exception):
    """This exception is used when a problem matching IMAGER data with Cloudsat / CALIPSO data has occured."""

    pass def get_matchups_from_data(cross, AM_PATHS, SETTINGS):
    """Find files and retrieve matchup from data."""

    # STEP 1 get imager files
    if SETTINGS['PPS_VALIDATION']:
        pps_files, imager_file, tobj = find_pps_cloud_files(cross, AM_PATHS, SETTINGS)
    if SETTINGS['CCI_CLOUD_VALIDATION']:
        imager_file, tobj = find_cci_cloud_file(cross, AM_PATHS)
    if SETTINGS['MAIA_VALIDATION']:
        imager_file, tobj = find_maia_cloud_file(cross, AM_PATHS)
    if SETTINGS['PATMOSX_VALIDATION']:
        imager_file, tobj = find_patmosx_cloud_file(cross, AM_PATHS)
    if SETTINGS['OCA_VALIDATION']:
        imager_file, tobj = find_oca_cloud_file(cross, AM_PATHS)
    if not imager_file:
        raise MatchupError("No imager file found!\ncross = " + str(cross))
    values = get_satid_datetime_orbit_from_fname(imager_file, SETTINGS, cross)
    date_time = values["date_time"]

    # Step 2 get truth satellite files
    truth_files = {}
    for truth in ['cloudsat', 'amsr', 'iss', 'synop', 'mora', 'cloudsat_lwp', 'calipso']:
        truth_files[truth] = None
        if (SETTINGS[truth.replace("_lwp", "").upper()+'_MATCHING'] and truth + '_file' in AM_PATHS.keys()):
            truth_files[truth] = find_truth_files(date_time, AM_PATHS, SETTINGS, values, truth=truth)
        elif not SETTINGS[truth.replace("_lwp", "").upper()+'_MATCHING']:
            logger.info("NO {truth} File, {truth} matching not requested "
                        "{truth}_MATCHING=False".format(truth=truth))
        elif truth + '_file' not in AM_PATHS.keys():
            logger.info("NO {truth}_file in atrain_match.cfg".format(truth=truth.lower()))
    # CALIPSO get some extra files:
    if truth_files['calipso'] is not None:
        extra_files = get_additional_calipso_files_if_requested(truth_files['calipso'], SETTINGS)
        calipso5km, calipso1km, calipso5km_aerosol = extra_files
    if (all(truth_files_i is None for truth_files_i in truth_files.values())):
        raise MatchupError(
            "Couldn't find any matching CALIPSO/CLoudSat/ISS data")

    # STEP 3 Read imager data:
    if (SETTINGS['PPS_VALIDATION']):
        cloudproducts = read_pps_data(pps_files, imager_file, SETTINGS)
        if os.path.isfile(SETTINGS['CNN_PCKL_PATH']):
            from atrain_match.utils.pps_prototyping_util import add_cnn_features_full
            cloudproducts.cnn_dict = add_cnn_features_full(cloudproducts.imager_channeldata,
                                                           cloudproducts,
                                                           SETTINGS)
    if (SETTINGS['CCI_CLOUD_VALIDATION']):
        cloudproducts = read_cloud_cci(imager_file)
        cloudproducts.satellite = values["satellite"]
    if (SETTINGS['MAIA_VALIDATION']):
        cloudproducts = read_cloud_maia(imager_file)
        cloudproducts.satellite = values["satellite"]
    if (SETTINGS['PATMOSX_VALIDATION']):
        cloudproducts = read_cloud_patmosx(imager_file, cross, SETTINGS)
        cloudproducts.satellite = values["satellite"]
    if (SETTINGS['OCA_VALIDATION']):
        cloudproducts = read_cloud_oca(imager_file)
        cloudproducts.satellite = values["satellite"]

    # STEP 4 get matchups
    # CloudSat
    cloudsat_matchup = None
    if ((SETTINGS['PPS_VALIDATION'] or SETTINGS['OCA_VALIDATION'])  and SETTINGS['CLOUDSAT_MATCHING'] and
            truth_files['cloudsat'] is not None):
        logger.info("Read CLOUDSAT data")
        cloudsat_matchup = get_cloudsat_matchups(truth_files['cloudsat'],
                                                 truth_files['cloudsat_lwp'],
                                                 cloudproducts, SETTINGS)
    # ISS:
    iss_matchup = None
    if (SETTINGS['PPS_VALIDATION'] and SETTINGS['ISS_MATCHING'] and
            truth_files['iss'] is not None):
        logger.info("Read ISS data")
        iss_matchup = get_iss_matchups(truth_files['iss'],
                                       cloudproducts, SETTINGS)
    # AMSR
    amsr_matchup = None
    if (SETTINGS['PPS_VALIDATION'] and SETTINGS['AMSR_MATCHING'] and
            truth_files['amsr'] is not None):
        logger.info("Read AMSR data")
        amsr_matchup = get_amsr_matchups(truth_files['amsr'],
                                         cloudproducts, SETTINGS)
    # SYNOP
    synop_matchup = None
    if (SETTINGS['PPS_VALIDATION'] and SETTINGS['SYNOP_MATCHING'] and
            truth_files['synop'] is not None):
        logger.info("Read SYNOP data")
        synop_matchup = get_synop_matchups(truth_files['synop'],
                                           cloudproducts, SETTINGS)
    # MORA
    mora_matchup = None
    if (SETTINGS['PPS_VALIDATION'] and SETTINGS['MORA_MATCHING'] and truth_files['mora'] is not None):
        logger.info("Read MORA data")
        mora_matchup = get_mora_matchups(truth_files['mora'],
                                         cloudproducts, SETTINGS)
    # CALIPSO:
    calipso_matchup = None
    if SETTINGS['CALIPSO_MATCHING'] and truth_files['calipso'] is not None:
        logger.info("Read CALIPSO data")
        calipso_matchup = get_calipso_matchups(truth_files['calipso'],
                                               values,
                                               cloudproducts,
                                               AM_PATHS, SETTINGS,
                                               calipso1km, calipso5km, calipso5km_aerosol)

    if calipso_matchup is None and SETTINGS['CALIPSO_REQUIRED']:
        raise MatchupError("No matches with CALIPSO.")
    elif cloudsat_matchup is None and SETTINGS['CLOUDSAT_REQUIRED']:
        raise MatchupError("No matches with CLOUSDAT.")
    elif iss_matchup is None and SETTINGS['ISS_REQUIRED']:
        raise MatchupError("No matches with ISS.")
    elif amsr_matchup is None and SETTINGS['AMSR_REQUIRED']:
        raise MatchupError("No matches with AMSR.")
    elif synop_matchup is None and SETTINGS['SYNOP_REQUIRED']:
        raise MatchupError("No matches with SYNOP.")
    elif mora_matchup is None and SETTINGS['MORA_REQUIRED']:
        raise MatchupError("No matches with MORA.")
    elif (calipso_matchup is None and cloudsat_matchup is None and iss_matchup is None
          and amsr_matchup is None and synop_matchup is None
          and mora_matchup is None):
        raise MatchupError("No matches with any truth.")

    # Get satellite name, time, and orbit number from imager_file
    date_time = values["date_time"]
    # basename = '_'.join(os.path.basename(imager_file).split('_')[:4])
    values = get_satid_datetime_orbit_from_fname(imager_file, SETTINGS, cross)
    basename = values["basename"]
    rematched_path = date_time.strftime(AM_PATHS['reshape_dir'].format(
        val_dir=config._validation_results_dir,
        satellite=values["satellite"],
        resolution=str(config.RESOLUTION),
    ))
    rematched_file = date_time.strftime(AM_PATHS['reshape_file'].format(
        satellite=values["satellite"],
        orbit=values["orbit"],
        resolution=str(config.RESOLUTION),
        instrument=INSTRUMENT.get(values["satellite"], 'avhrr'),
        atrain_datatype="atrain_datatype",
        extrai=values.get("extrai", "")
    ))
    rematched_file_base = rematched_path + rematched_file

    # Create directories if they don't exist yet
    if not os.path.exists(os.path.dirname(rematched_path)):
        logger.info("Creating dir %s:", rematched_path)
        os.makedirs(os.path.dirname(rematched_path))

    for matchup, name in zip([cloudsat_matchup, iss_matchup, amsr_matchup,
                              synop_matchup, mora_matchup, calipso_matchup],
                             ['CloudSat', 'ISS', 'AMSR-E',
                              'SYNOP', 'MORA', 'CALIPSO']):
        if matchup is None:
            continue
        # add modis lvl2
        if SETTINGS['MATCH_MODIS_LVL2']:
            from atrain_match.cloudproducts.read_modis_products import add_modis_06
            if matchup.imager_instrument in ['modis']:
                matchup = add_modis_06(matchup, AM_PATHS, cross)
        if SETTINGS['ADD_NWP']:
            import pps_nwp
            from atrain_match.libs.extract_imager_along_track import _interpolate_height_and_temperature_from_pressure
            nwp_file = find_closest_nwp_file(cloudproducts, AM_PATHS,
                                             values, SETTINGS)
            logger.debug(nwp_file)
            if pps_nwp is None:
                continue
            gribfile = pps_nwp.GRIBFile(nwp_file, (matchup.imager.longitude,
                                                   matchup.imager.latitude))
            setattr(matchup.imager, "nwp_height",
                    gribfile.get_gh_vertical()[0, :, :].astype(np.float32).transpose())
            setattr(matchup.imager, "nwp_surface_h",
                    gribfile.get_gh_surface()[:].astype(np.float32))
            setattr(matchup.imager, "nwp_temperature",
                    gribfile.get_t_vertical()[0, :, :].astype(np.float32).transpose())
            setattr(matchup.imager, "nwp_h2m",
                    gribfile.get_h_2meter()[:].astype(np.float32))
            setattr(matchup.imager, "nwp_t2m",
                    gribfile.get_t_2meter()[:].astype(np.float32))
            setattr(matchup.imager, "nwp_u10m",
                    gribfile.get_u_10meter()[:].astype(np.float32))
            setattr(matchup.imager, "nwp_v10m",
                    gribfile.get_v_10meter()[:].astype(np.float32))
            # get pressure variables in hPs
            field = gribfile.get_p_vertical()
            if field.units == 'Pa':
                field = 0.01*field[:]
                field.units = 'hPa'
            matchup.imager.nwp_pressure = field[0, :, :].astype(np.float32).transpose()
            field = gribfile.get_p_surface()
            if field.units == 'Pa':
                field = 0.01*field[:]
                field.units = 'hPa'
            matchup.imager.nwp_psur = field[:].astype(np.float32).transpose()
            data = _interpolate_height_and_temperature_from_pressure(matchup.imager, 440)
            setattr(matchup.imager, 'nwp_h440', data)
            data = _interpolate_height_and_temperature_from_pressure(matchup.imager, 680)
            setattr(matchup.imager, 'nwp_h680', data)
            if SETTINGS['OCA_VALIDATION'] and matchup.imager.ctth_height is None:
                data = _interpolate_height_and_temperature_from_pressure(matchup.imager, None,
                                                                         list_of_levels=matchup.imager.ctth_pressure)
                data[matchup.imager.ctth_pressure<0] = -9
                setattr(matchup.imager, 'ctth_height', data)

    # add additional vars to cloudsat and calipso objects and print them to file:
    cloudsat_matchup, calipso_matchup = add_additional_clousat_calipso_index_vars(cloudsat_matchup, calipso_matchup)
    cloudsat_matchup, calipso_matchup, iss_matchup = add_elevation_corrected_imager_ctth(
        cloudsat_matchup, calipso_matchup, iss_matchup, SETTINGS)

    # imager_name
    imager_obj_name = 'pps'
    if SETTINGS['CCI_CLOUD_VALIDATION']:
        imager_obj_name = 'cci'
    if SETTINGS['OCA_VALIDATION']:
        imager_obj_name = 'oca'
    if SETTINGS['MAIA_VALIDATION']:
        imager_obj_name = 'maia'
    if SETTINGS['PATMOSX_VALIDATION']:
        imager_obj_name = 'patmosx'

    # write matchups
    for matchup, name in zip([cloudsat_matchup, iss_matchup, amsr_matchup,
                              synop_matchup, mora_matchup, calipso_matchup],
                             ['CloudSat', 'ISS', 'AMSR-E',
                              'SYNOP', 'MORA', 'CALIPSO']):
        if matchup is None:
            logger.debug("No {:s} Match File created".format(name))
        else:
            truth_sat = matchup.truth_sat
            match_file = rematched_file_base.replace(
                'atrain_datatype', truth_sat)
            write_truth_imager_match_obj(match_file, matchup,
                                         SETTINGS,
                                         imager_obj_name=imager_obj_name)

    # no longer return the matchup data?
    return {'cloudsat': cloudsat_matchup,
            'calipso': calipso_matchup,
            'iss': iss_matchup,
            'amsr': amsr_matchup,
            'synop': synop_matchup,
            'mora': mora_matchup,
            'basename': basename,
            'values': values} def find_main_cloudproduct_file(cross, filedir_pattern, filename_pattern, values=None):
    """Find main cloudproduct file with lat/lon (inner function)."""
    if values is None:
        values = {}

    # tlist = get_time_list(cross.time,
    #                      (timedelta(seconds=0), timedelta(seconds=0)),
    #                      delta_t_in_seconds=60)
    tlist = [cross.time]
    values["satellite"] = cross.satellite1
    found_dir = None
    no_files_found = True
    checked_dir = {}
    for tobj in tlist:
        # print values
        imager_dir = insert_info_in_filename_or_path(filedir_pattern,
                                                     values, datetime_obj=tobj)
        # print imager_dir
        if os.path.exists(imager_dir):
            found_dir = imager_dir
            if imager_dir not in checked_dir.keys():
                checked_dir[imager_dir] = 1
                logger.info("Found directory %s", found_dir)
        else:
            if not found_dir and imager_dir not in checked_dir.keys():
                checked_dir[imager_dir] = 1
                logger.info("This directory does not exist, pattern: %s", imager_dir)
            continue
        file_pattern = insert_info_in_filename_or_path(filename_pattern,
                                                       values, datetime_obj=tobj)
        files = glob(os.path.join(found_dir, file_pattern))
        if len(files) > 0:
            no_files_found = False
            logger.info("Found file: %s", os.path.basename(str(files[0])))
            return files[0], tobj
    if not found_dir:
        return None, None
    if no_files_found:
        logger.info("Found no files for patterns of type: %s", file_pattern)
    return None, None def insert_info_in_filename_or_path(file_or_name_path, values, datetime_obj=None):
    """Fill placeholders in filename/pathnames from atrain_match.cfg with correct data."""
    # file_or_name_path = file_or_name_path.format(**values)
    satellite = values.get("satellite", "*")
    file_or_name_path = file_or_name_path.format(
        ctth_type=values.get("ctth_type", ""),
        satellite=satellite,
        orbit=values.get("orbit", "*"),
        instrument=INSTRUMENT.get(satellite, "avhrr"),
        resolution=config.RESOLUTION,
        lines_lines=values.get("lines_lines", "*"),
        val_dir=config._validation_results_dir,
        extrai=values.get('extrai', ""),
        plus_hours=values.get("plus_hours", ""),
        year=values.get('year', "unknown"),
        month=values.get('month', "unknown"),
        mode=values.get('mode', "unknown"),
        # min_opt_depth=values.get('min_opt_depth', ""),
        atrain_datatype=values.get("atrain_datatype", "atrain_datatype"))
    if datetime_obj is None:
        return file_or_name_path
    name = datetime_obj.strftime(file_or_name_path)
    return name def add_additional_clousat_calipso_index_vars(match_clsat, match_calipso):
    # add cloudsat modisflag to calipso obj

    if match_clsat is not None and match_calipso is not None:
        match_calipso.calipso.cal_modis_cflag = None
        # map cloudsat to calipso and the other way around!
        from atrain_match.utils.match import match_lonlat
        source = (match_clsat.cloudsat.longitude.astype(np.float64).reshape(-1, 1),
                  match_clsat.cloudsat.latitude.astype(np.float64).reshape(-1, 1))
        target = (match_calipso.calipso.longitude.astype(np.float64).reshape(-1, 1),
                  match_calipso.calipso.latitude.astype(np.float64).reshape(-1, 1))
        mapper, dummy = match_lonlat(source, target, radius_of_influence=1000, n_neighbours=1)
        match_calipso.calipso.cloudsat_index = mapper.rows.filled(config.NODATA).ravel()
        target = (match_clsat.cloudsat.longitude.astype(np.float64).reshape(-1, 1),
                  match_clsat.cloudsat.latitude.astype(np.float64).reshape(-1, 1))
        source = (match_calipso.calipso.longitude.astype(np.float64).reshape(-1, 1),
                  match_calipso.calipso.latitude.astype(np.float64).reshape(-1, 1))
        mapper, dummy = match_lonlat(source, target, radius_of_influence=1000, n_neighbours=1)
        match_clsat.cloudsat.calipso_index = mapper.rows.filled(config.NODATA).ravel()

        # Transfer CloudSat MODIS cloud flag to CALIPSO representation
        index = match_calipso.calipso.cloudsat_index.copy()
        index[index < 0] = 0
        match_calipso.calipso.cal_modis_cflag = np.where(
            match_calipso.calipso.cloudsat_index >= 0,
            match_clsat.cloudsat.MODIS_cloud_flag[index],
            -9)

    if match_clsat is not None and match_calipso is not None:
        index = match_clsat.cloudsat.calipso_index.copy()
        index[index < 0] = 0
        for var_2d_name in ['feature_classification_flags',
                            'layer_base_altitude',
                            'layer_top_altitude',
                            'feature_optical_depth_532',
                            'feature_optical_depth_532_5km']:
            if hasattr(match_calipso.calipso, var_2d_name):
                data_calipso = getattr(match_calipso.calipso, var_2d_name)
                if data_calipso is None:
                    continue
                temp_data = np.array([np.where(
                    match_clsat.cloudsat.calipso_index >= 0,
                    data_calipso[index, i], -9) for i in range(data_calipso.shape[1])])
                setattr(match_clsat.cloudsat, 'calipso_{:s}'.format(var_2d_name), temp_data.transpose())

        for var_1d_name in ['column_optical_depth_tropospheric_aerosols_532',
                            'column_optical_depth_tropospheric_aerosols_532_5km',
                            'column_optical_depth_aerosols_532',
                            'column_optical_depth_aerosols_532_5km']:

            if hasattr(match_calipso.calipso, var_1d_name):
                data_calipso = getattr(match_calipso.calipso, var_1d_name)
                if data_calipso is None:
                    continue
                temp_data = np.where(match_clsat.cloudsat.calipso_index >= 0, data_calipso[index], -9)
                setattr(match_clsat.cloudsat, 'calipso_{:s}'.format(var_1d_name), temp_data)

        """
        match_clsat.cloudsat.calipso_feature_classification_flags= np.where(
            match_clsat.cloudsat.calipso_index >= 0,
            match_calipso.calipso.feature_classification_flags[index, 0],
            -9)
        # first bas layer use height not pressure as cloudsat uses height
        match_clsat.cloudsat.calipso_layer_base_altitude = np.where(
            match_clsat.cloudsat.calipso_index >= 0,
            match_calipso.calipso.layer_base_altitude[index, 0],
            -9)
        for layer in range(1, 10):
            match_clsat.cloudsat.calipso_layer_base_altitude = np.where(
                np.logical_and(np.logical_and(match_clsat.cloudsat.calipso_index >= 0,
                                              match_calipso.calipso.layer_base_altitude[index, layer]>0),
                               match_calipso.calipso.layer_base_altitude[index, layer]<
                               match_calipso.calipso.layer_base_altitude[index, layer-1]),

                match_calipso.calipso.layer_base_altitude[index, layer],
            match_clsat.cloudsat.calipso_layer_base_altitude)
        match_clsat.cloudsat.calipso_layer_base_altitude[match_clsat.cloudsat.calipso_layer_base_altitude<-999] = -9.0
        # first bas layer use height not pressure as cloudsat uses height
        match_clsat.cloudsat.calipso_layer_top_altitude = np.where(
            match_clsat.cloudsat.calipso_index >= 0,
            match_calipso.calipso.layer_top_altitude[index, 0],
            -9)
        for layer in range(1, 10):
            match_clsat.cloudsat.calipso_layer_top_altitude = np.where(
                np.logical_and(np.logical_and(match_clsat.cloudsat.calipso_index >= 0,
                                              match_calipso.calipso.layer_top_altitude[index, layer]>0),
                               match_calipso.calipso.layer_top_altitude[index, layer]>
                               match_calipso.calipso.layer_top_altitude[index, layer-1]),

                match_calipso.calipso.layer_top_altitude[index, layer],
            match_clsat.cloudsat.calipso_layer_top_altitude)
        match_clsat.cloudsat.calipso_layer_top_altitude[match_clsat.cloudsat.calipso_layer_top_altitude<-999] = -9.0
    """
    return match_clsat, match_calipso def add_elevation_corrected_imager_ctth(match_clsat, match_calipso, match_iss, SETTINGS):
    # # Cloudsat # #
    if match_clsat is None or match_clsat.imager.ctth_height is None:
        pass
    elif match_clsat.imager.imager_ctth_m_above_seasurface is None:
        # First make sure that PPS cloud top heights are converted to height
        # above sea level just as CloudSat height are defined. Use
        # corresponding DEM data.
        elevation = np.where(np.less_equal(match_clsat.cloudsat.elevation, 0),
                             0, match_clsat.cloudsat.elevation)
        num_csat_data_ok = len(match_clsat.cloudsat.elevation)
        logger.debug("Length of CLOUDSAT array: %d", num_csat_data_ok)
        imager_ctth_m_above_seasurface = match_clsat.imager.ctth_height.copy()
        # import pdb;pdb.set_trace()
        if SETTINGS["CCI_CLOUD_VALIDATION"] or SETTINGS["PATMOSX_VALIDATION"]:
            # ctth already relative mean sea level
            imager_ctth_m_above_seasurface = match_calipso.imager.ctth_height
        else:  # ctth relative topography
            got_height = imager_ctth_m_above_seasurface >= 0
            imager_ctth_m_above_seasurface[got_height] += elevation[got_height]*1.0
        match_clsat.imager.imager_ctth_m_above_seasurface = imager_ctth_m_above_seasurface
        if num_csat_data_ok == 0:
            logger.info("Processing stopped: Zero lenght of matching arrays!")
            raise ProcessingError("Zero lenght of matching arrays!")
    # # Calipso # #
    # First make sure that PPS cloud top heights are converted to height above sea level
    # just as CALIPSO heights are defined. Use corresponding DEM data.
    if match_calipso is None or match_calipso.imager.ctth_height is None:
        pass
    elif match_calipso.imager.imager_ctth_m_above_seasurface is None:
        cal_elevation = np.where(np.less_equal(match_calipso.calipso.elevation, 0),
                                 0, match_calipso.calipso.elevation)
        num_cal_data_ok = len(match_calipso.calipso.elevation)
        logger.debug("Length of CALIOP array: %d", num_cal_data_ok)
        imager_ctth_m_above_seasurface = match_calipso.imager.ctth_height.copy()
        logger.debug("CCI_CLOUD_VALIDATION %s", str(SETTINGS["CCI_CLOUD_VALIDATION"]))
        if SETTINGS["CCI_CLOUD_VALIDATION"]:
            # ctth relative mean sea level
            imager_ctth_m_above_seasurface = match_calipso.imager.ctth_height
        else:  # ctth relative topography
            got_height = imager_ctth_m_above_seasurface >= 0
            imager_ctth_m_above_seasurface[got_height] += cal_elevation[got_height]*1.0
        match_calipso.imager.imager_ctth_m_above_seasurface = imager_ctth_m_above_seasurface
    if match_iss is None or match_iss.imager.ctth_height is None:
        pass
    elif match_iss.imager.imager_ctth_m_above_seasurface is None:
        iss_elevation = np.where(np.less_equal(match_iss.iss.elevation, 0),
                                 0, match_iss.iss.elevation)
        num_iss_data_ok = len(match_iss.iss.elevation)
        logger.info("Length of ISS array: %d", num_iss_data_ok)
        imager_ctth_m_above_seasurface = match_iss.imager.ctth_height.copy()
        if SETTINGS["CCI_CLOUD_VALIDATION"]:
            # ctth relative mean sea level
            pass
        else:  # ctth relative topography
            got_height = imager_ctth_m_above_seasurface >= 0
            imager_ctth_m_above_seasurface[got_height] += iss_elevation[got_height]*1.0
        match_iss.imager.imager_ctth_m_above_seasurface = imager_ctth_m_above_seasurface
    return match_clsat, match_calipso, match_iss def add_validation_ctth_cloudsat(cloudsat):
    # CLOUDSAT VALIDATION HEIGHT!
    # The restriction to use only pixels where imager and cloudsat
    # both is cloudy is done later. This makes it possbile to find also
    # POD-cloudy for different cloud heights
    LARGE_POSITIVE = 99999.0
    validation_height = -9 + 0 * np.zeros(cloudsat.latitude.shape)
    validation_height_base = LARGE_POSITIVE + 0 * np.zeros(cloudsat.latitude.shape)
    for i in range(125):
        height = cloudsat.Height[:, i]
        cmask_ok = cloudsat.CPR_Cloud_mask[:, i]
        top_height = height + 120
        # top_height[height<240*4] = -9999 # Do not use not sure why these were not used Nina 20170317
        is_cloudy = cmask_ok > config.CLOUDSAT_CLOUDY_THR
        top_height[~is_cloudy] = -9999
        validation_height[validation_height < top_height] = top_height[
            validation_height < top_height]
        cloudsat.validation_height = validation_height
        update_base = np.logical_and(top_height > 0, validation_height_base > top_height)
        validation_height_base[update_base] = top_height[update_base]
        cloudsat.validation_height_base = validation_height_base
    cloudsat.validation_height_base[cloudsat.validation_height_base >= LARGE_POSITIVE] = -9
    return cloudsat def add_cloudsat_cloud_fraction(cloudsat):
    cloudsat_cloud_mask = cloudsat.CPR_Cloud_mask
    cloudsat_cloud_mask = np.greater_equal(cloudsat_cloud_mask,
                                           config.CLOUDSAT_CLOUDY_THR)
    cloudsat_cloud_fraction = np.zeros(cloudsat.latitude.shape[0])
    sum_cloudsat_cloud_mask = np.sum(cloudsat_cloud_mask, axis=1)
    if len(sum_cloudsat_cloud_mask) != (len(cloudsat_cloud_fraction)):
        raise ValueError('Boolean index-array should have same lenght as array!')
    cloudsat_cloud_fraction[sum_cloudsat_cloud_mask > 2] = 1  # requires at least two cloudy bins
    cloudsat.cloud_fraction = cloudsat_cloud_fraction
    return cloudsat def optical_depth_height_filtering(calipso, min_optical_depth, use_old_method=False,
                                   limit_ctop=0.2):
    new_cloud_top = np.ones(calipso.layer_top_altitude.shape)*config.NODATA*1.0
    new_cloud_base = np.ones(calipso.layer_base_altitude.shape)*config.NODATA*1.0
    new_cloud_top_pressure = np.ones(calipso.layer_top_altitude.shape)*config.NODATA*1.0
    new_cloud_base_pressure = np.ones(calipso.layer_base_altitude.shape)*config.NODATA*1.0
    distance_down_in_cloud_we_see = np.ones(calipso.layer_base_altitude.shape)*config.NODATA*1.0
    new_cloud_fraction = np.zeros(calipso.cloud_fraction.shape)
    new_fcf = np.ones(calipso.feature_classification_flags.shape).astype(calipso.feature_classification_flags.dtype)
    depthsum = 0.0*new_cloud_fraction.copy()
    already_detected = depthsum > 99999  # None from the start! All False
    already_detected[calipso.layer_top_altitude[:, 0] < 0] = True  # don't bother with clear pixels
    # to_thin_to_see_at_all = calipso.total_optical_depth_5km < min_optical_depth
    N10 = new_cloud_top.shape[1]
    for layer_j in range(N10):
        # Filter OPTICAL_LIMIT_CLOUD_TOP down in EACH layer
        # Notice even if top of layer i is always above layer i-1,
        # it is NOT the case that base of layer i is always above layer i-1!
        # We could have the situation:
        # layer 0 base at 2km top at 0km optical thickness 0.1
        # layer 1 base at 8km top at 9km optical thickness 10
        # Or the situation:
        # layer 0 base at 2km top at 0km optical thickness 0.5
        # layer 1 base at 8km top at 9km optical thickness 0.5
        this_layer_optical_depth = calipso.feature_optical_depth_532[:, layer_j].ravel()
        cloud_base = calipso.layer_base_altitude[:, layer_j]
        cloud_top = calipso.layer_top_altitude[:, layer_j]
        geometrical_distance_cloud = cloud_top-cloud_base
        geometrical_distance_cloud[cloud_base < 0] = 0
        geometrical_distance_cloud[cloud_top < 0] = 0
        fraction_into_cloud = (limit_ctop*1.0)/this_layer_optical_depth
        fraction_into_cloud[fraction_into_cloud < 0] = 0
        fraction_into_cloud[fraction_into_cloud > 1.0] = 1.0
        if use_old_method:
            fraction_into_cloud = 0.5
        # For KGs method set fraction_into_cloud = 0.5 always
        distance_down_in_cloud_we_see[:, layer_j] = geometrical_distance_cloud*fraction_into_cloud
    for layer_j in range(N10):
        # total_optical_depth_layers_above = depthsum.copy()
        this_layer_optical_depth = calipso.feature_optical_depth_532[:, layer_j].ravel()
        depthsum += this_layer_optical_depth
        # this_is_the_top_most_layer_we_detect
        update = np.logical_and(
            depthsum >= min_optical_depth,
            np.equal(already_detected, False))
        new_cloud_top[update, 0:(N10-layer_j)] = (calipso.layer_top_altitude[update, layer_j:] -
                                                  distance_down_in_cloud_we_see[update, layer_j:])
        new_cloud_base[update, 0:(N10-layer_j)] = calipso.layer_base_altitude[update, layer_j:]
        new_cloud_top_pressure[update, 0:(N10-layer_j)] = (calipso.layer_top_pressure[update, layer_j:])
        new_cloud_base_pressure[update, 0:(N10-layer_j)] = calipso.layer_base_pressure[update, layer_j:]
        new_fcf[update, 0:(N10-layer_j)] = calipso.feature_classification_flags[update, layer_j:]
        new_cloud_fraction[update] = calipso.cloud_fraction[update]
        already_detected[update] = True
    if use_old_method:
        pass
    else:
        for layer_j in range(N10):
            filtered_to_low = new_cloud_top[:, 0] < new_cloud_top[:, layer_j]
            new_cloud_top[filtered_to_low, 0] = new_cloud_top[filtered_to_low, layer_j]

    lowest_positive_cloud_base = calipso.layer_base_altitude.copy()
    lowest_positive_cloud_base[lowest_positive_cloud_base < 0] = 99999
    lowest_positive_cloud_base = np.min(lowest_positive_cloud_base, axis=1)
    lowest_positive_cloud_base[lowest_positive_cloud_base == 99999] = -9

    new_validation_height = new_cloud_top[:, 0].copy()
    new_validation_height[new_validation_height >= 0] = new_validation_height[new_validation_height >= 0]*1000
    new_validation_height[new_validation_height < 0] = -9

    detection_height = new_validation_height.copy()
    detection_height[detection_height < 0] = lowest_positive_cloud_base[detection_height < 0] * 1000.0  #ok in val h?
    detection_height[detection_height < 0] = -9
    return (new_cloud_top, new_cloud_base, new_cloud_fraction, new_fcf, new_validation_height,
            new_cloud_top_pressure, new_cloud_base_pressure, detection_height) def check_total_optical_depth_and_warn(match_calipso):
    obj = match_calipso.calipso
    if (obj.total_optical_depth_5km is not None and
            (obj.total_optical_depth_5km < obj.feature_optical_depth_532_top_layer_5km).any()):
        print("warning {:d}".format(len(obj.total_optical_depth_5km)))
        # print len(obj.total_optical_depth_5km[badPix])
        # print obj.total_optical_depth_5km[badPix]
        # print obj.feature_optical_depth_532_top_layer_5km[badPix]
        # print diff[badPix]
        # print obj.number_layers_found[badPix]
        # if obj.detection_height_5km is not None:
        #    print obj.detection_height_5km[badPix]
        # print np.where(badPix)
        # print obj.layer_top_altitude[0, badPix]
        # print obj.layer_base_altitude[0, badPix] def add_validation_ctth_calipso(calipso):
    calipso.validation_height = calipso.layer_top_altitude[:, 0].copy()
    calipso.validation_height[calipso.validation_height >= 0] = (
        1000.0 * calipso.validation_height[calipso.validation_height >= 0])
    calipso.validation_height[calipso.validation_height < 0] = -9
    return calipso def detection_height_filtering(match_calipso):
    if match_calipso.calipso.detection_height_5km is None:
        logger.warning("Reprocess matchups with CALCULATE_DETECTION_HEIGHT_FROM_5KM_DATA = True")
        raise ProcessingError("No detection_height_5km in rehsaped file!")

    # Should all layers be updated???
    new_cloud_tops = np.where(
        match_calipso.calipso.layer_top_altitude[:, 0]*1000 > match_calipso.calipso.detection_height_5km,
        match_calipso.calipso.detection_height_5km,
        match_calipso.calipso.layer_top_altitude[:, 0]*1000)
    new_cloud_tops = np.where(
        match_calipso.calipso.layer_base_altitude[:, 0]*1000 > match_calipso.calipso.detection_height_5km,
        match_calipso.calipso.layer_base_altitude[:, 0]*1000,
        new_cloud_tops)
    clouds_to_update = np.logical_and(
        match_calipso.calipso.layer_top_altitude[:, 0]*1000 > match_calipso.calipso.detection_height_5km,
        np.not_equal(match_calipso.calipso.detection_height_5km, -9))  # Dont update clouds only in 1km
    return np.where(clouds_to_update,
                    new_cloud_tops,
                    match_calipso.calipso.validation_height) def set_thin_to_clear_filtering_1km(match_calipso, SETTINGS):
    isThin_clouds = np.logical_and(
        match_calipso.calipso.total_optical_depth_5km < SETTINGS['OPTICAL_DETECTION_LIMIT'],
        np.not_equal(match_calipso.calipso.total_optical_depth_5km, -9))
    # >0.0 important. Some clouds are missing in 5km data set but present in 1km data set!
    set_to_clear = np.logical_and(
        match_calipso.calipso.number_layers_found > 0,
        isThin_clouds)
    cloud_fraction = match_calipso.calipso.cloud_fraction
    validation_height = match_calipso.calipso.validation_height
    cloud_fraction[set_to_clear] = 0.00001
    validation_height[set_to_clear] = -9
    return cloud_fraction, validation_height def calculate_statistics(mode, statfilename, match_calipso, match_clsat, match_iss, match_amsr, match_synop,
                         SETTINGS,
                         dnt_flag=None):
    """Calculate all requested statistics, for all matches with one imager cloudproduct file (main function)."""

    def get_day_night_subset(match_obj, val_subset, SETTINGS):
        (no_qflag, night_flag, twilight_flag,
         day_flag, all_dnt_flag) = get_day_night_info(match_obj, SETTINGS)

        if dnt_flag is None:
            logger.debug('dnt_flag = %s', 'ALL PIXELS')
            dnt_subset = np.logical_and(val_subset, all_dnt_flag)
        elif dnt_flag.upper() == 'DAY':
            logger.debug('dnt_flag = %s', dnt_flag.upper())
            dnt_subset = np.logical_and(val_subset, day_flag)
        elif dnt_flag.upper() == 'NIGHT':
            logger.debug('dnt_flag = %s', dnt_flag.upper())
            dnt_subset = np.logical_and(val_subset, night_flag)
        elif dnt_flag.upper() == 'TWILIGHT':
            logger.debug('dnt_flag = %s', dnt_flag.upper())
            dnt_subset = np.logical_and(val_subset, twilight_flag)
        else:
            raise ProcessingError("Unknown DNT-flag %s" % (dnt_flag.upper()))
        return dnt_subset

    if match_clsat is not None:
        logger.info("Cloudsat Statistics")
        val_subset = get_subset_for_mode(match_clsat, mode)
        if val_subset is not None:
            val_subset = get_day_night_subset(match_clsat, val_subset, SETTINGS)
            statfile = open(statfilename.replace('xxx', 'cloudsat'), "w")
            if match_clsat.cloudsat.all_arrays['cloud_fraction'] is not None:
                low_medium_high_class = get_cloudsat_low_medium_high_classification(match_clsat)
                print_main_stats(match_clsat, statfile)
                print_cmask_stats(match_clsat, statfile, val_subset, SETTINGS)
                print_cmask_prob_stats(match_clsat, statfile, val_subset, SETTINGS)
                print_modis_stats(match_clsat, statfile, val_subset, match_clsat.cloudsat.MODIS_cloud_flag, SETTINGS)
                print_stats_ctop(match_clsat, statfile, val_subset, low_medium_high_class, SETTINGS)
            if match_clsat.cloudsat.all_arrays['liq_water_path_gm2'] is not None:
                print_cpp_lwp_stats(match_clsat, statfile, val_subset)
            statfile.close()

    if match_calipso is not None:
        logger.info("Calipo Statistics")
        val_subset = get_subset_for_mode(match_calipso, mode)
        if val_subset is not None:
            statfile = open(statfilename.replace('xxx', 'calipso'), "w")
            low_medium_high_class = get_calipso_low_medium_high_classification(match_calipso)
            # semi_flag, opaque_flag = get_semi_opaque_info(match_calipso)
            val_subset = get_day_night_subset(match_calipso, val_subset, SETTINGS)
            print_main_stats(match_calipso, statfile)
            print_cmask_stats(match_calipso, statfile, val_subset, SETTINGS)
            print_cmask_prob_stats(match_calipso, statfile, val_subset, SETTINGS)
            print_modis_stats(match_calipso, statfile, val_subset, match_calipso.calipso.cal_modis_cflag, SETTINGS)
            print_calipso_stats_ctype(match_calipso, statfile, val_subset, low_medium_high_class, SETTINGS)
            print_stats_ctop(match_calipso, statfile, val_subset, low_medium_high_class, SETTINGS)
            print_cpp_stats(match_calipso, statfile, val_subset, SETTINGS)
            statfile.close()

    if match_iss is not None:
        val_subset = get_subset_for_mode(match_iss, mode)
        if val_subset is not None:
            statfile = open(statfilename.replace('xxx', 'iss'), "w")
            val_subset = get_day_night_subset(match_iss, val_subset, SETTINGS)
            print_main_stats(match_iss, statfile)
            print_cmask_stats(match_iss, statfile, val_subset, SETTINGS)
            print_cmask_prob_stats(match_iss, statfile, val_subset, SETTINGS)
            # print_calipso_stats_ctype(match_iss, statfile, val_subset, cal_vert_feature)
            print_stats_ctop(match_iss, statfile, val_subset, None, SETTINGS)
            statfile.close()

    if match_amsr is not None:
        logger.info("AMSR-E Statistics")
        val_subset = get_subset_for_mode(match_amsr, mode)
        if val_subset is not None:
            val_subset = get_day_night_subset(match_amsr, val_subset, SETTINGS)
            statfile = open(statfilename.replace('xxx', 'amsr'), "w")
            print_main_stats(match_amsr, statfile)
            print_cpp_lwp_stats(match_amsr, statfile, val_subset)
            statfile.close()

    if match_synop is not None:
        logger.info("SYNOP Statistics")
        val_subset = get_subset_for_mode(match_synop, mode)
        if val_subset is not None:
            val_subset = get_day_night_subset(match_synop, val_subset, SETTINGS)
            statfile = open(statfilename.replace('xxx', 'synop'), "w")
            print_main_stats(match_synop, statfile)
            print_cmask_stats(match_synop, statfile, val_subset, SETTINGS)
            print_cmask_prob_stats(match_synop, statfile, val_subset, SETTINGS)
            statfile.close() def plot_satellite_trajectory(longitude,
                              latitude,
                              trajectoryname,
                              area_config_file,
                              fig_type='eps',
                              **options):
    """Plot a trajectory of geolocations (lon, lat) on a map."""
    import numpy as np

    if 'trajectory_plot_area' in options:
        area_id = options['trajectory_plot_area']
    else:
        area_id = "pc_world"  # Global area
    print(area_id)

    track = np.ones(latitude.shape)

    import pyresample as pr
    swath_def = pr.geometry.SwathDefinition(longitude, latitude)

    area_def = pr.utils.load_area(area_config_file, area_id)

    result = pr.kd_tree.resample_nearest(swath_def, track, area_def,
                                         radius_of_influence=50 * 1000, fill_value=None)

    import matplotlib.pyplot as plt
    import pylab
    fig = plt.figure()
    ax = fig.add_subplot(111)

    bmap = pr.plot.area_def2basemap(area_def)
    bmng = bmap.bluemarble()
    # col = bmap.imshow(result, cmap=pylab.get_cmap('rainbow'), origin='upper')
    col = bmap.imshow(result, origin='upper', cmap=pylab.get_cmap('autumn'))
    # plt.title('Calipso matchup track')

    for figt in fig_type:
        figname = trajectoryname + '.' + figt
        plt.savefig(figname, bbox_inches='tight')

    del plt
    return def plot_cal_clsat_imager_time_diff(match_clsat,
                                    match_calipso,
                                    plotpath, basename,
                                    resolution, file_type='png',
                                    **options):
    """Plot time difference between imager and CALIOP, CPR (CloudSat)."""
    if 'instrument' in options:
        instrument = options['instrument']
    else:
        instrument = 'imager'
    pixel_position = np.arange(match_calipso.calipso.latitude.shape[0])
    cal_diff_sec_1970 = match_calipso.diff_sec_1970 / 60.0
    # Plot time diff
    fig = plt.figure()
    ax = fig.add_subplot(111)
    maxvalue = np.nanmax(cal_diff_sec_1970)
    minvalue = np.nanmin(cal_diff_sec_1970)
    title = "Time Difference CALIPSO - %s" % instrument.upper()
    ylabel_str = "Time diff (CALIPSO - %s)[min]" % instrument.upper()
    if match_clsat is not None:
        clsat_diff_sec_1970 = match_clsat.diff_sec_1970 / 60.0
        title = "Time Difference Between %s and CloudSat/CALIPSO" % instrument.upper()
        ylabel_str = "Time diff (CALIPSO/CloudSat - %s)[min]" % instrument.upper()
        maxvalue = np.max([np.nanmax(clsat_diff_sec_1970), maxvalue])
        minvalue = np.min([np.nanmin(clsat_diff_sec_1970), minvalue])

        biggest_Cloudsat_diff = np.nanmax(np.abs(clsat_diff_sec_1970 ))
        # For 5km data many cloudsat index are -9 as cloudsat has smalles res
        # For now just exclude those
        use = match_clsat.cloudsat.calipso_index > 0
        ax.plot(pixel_position[match_clsat.cloudsat.calipso_index[use]],
                clsat_diff_sec_1970[use], 'r+',
                label="CloudSat (max time diff %.2f min)" % (biggest_Cloudsat_diff))

    biggest_Calipso_diff = np.nanmax(np.abs(cal_diff_sec_1970 ))
    ax.set_title(title)
    ax.set_xlabel("Track Position")
    ax.set_ylabel(ylabel_str)
    ax.set_ylim(minvalue - 1, maxvalue + 1)
    ax.plot(pixel_position, cal_diff_sec_1970, "g",
            label="CALIPSO (max time diff %.2f min)" % (biggest_Calipso_diff))
    ax.legend(numpoints=4)

    if isinstance(file_type, str):
        fig.savefig("%s/%skm_%s_time_diff.%s" % (plotpath,
                                                 resolution, basename, file_type))
    else:
        for filetype in file_type:
            fig.savefig("%s/%skm_%s_time_diff.%s" % (plotpath,
                                                     resolution, basename, filetype))
    return def plot_cal_clsat_geoprof_imager(match_clsat,
                                  match_calipso,
                                  imager_ctth_m_above_seasurface,
                                  plotpath,
                                  basename,
                                  mode,
                                  file_type='png',
                                  **options):
    """Plot imager co-located data on CALIOP and CPR (Cloudsat) track."""
    instrument = 'imager'
    if 'instrument' in options:
        instrument = options['instrument']
    MAXHEIGHT = None
    if 'MAXHEIGHT' in options:
        MAXHEIGHT = options["MAXHEIGHT"]
    caliop_height = match_calipso.calipso.layer_top_altitude * 1000
    caliop_base = match_calipso.calipso.layer_base_altitude * 1000
    calipso_val_h = match_calipso.calipso.validation_height
    caliop_base[caliop_base < 0] = -9
    caliop_height[caliop_height < 0] = -9
    if RESOLUTION == 5 and match_clsat is not None:
        # Cloudsat has highest res
        pixel_position = np.arange(match_clsat.cloudsat.latitude.shape[0])[match_calipso.calipso.cloudsat_index]
    else:    
        pixel_position = np.arange(match_calipso.calipso.latitude.shape[0])    
    # Calculates Hihest Cloud Top
    if MAXHEIGHT is None:
        maxheight_calipso = np.nanmax(caliop_height)
        maxheight_imager = np.nanmax(imager_ctth_m_above_seasurface)
        max_height_sat = np.max([maxheight_calipso, maxheight_imager])
        maxheight = max_height_sat + 1000
    else:
        maxheight = MAXHEIGHT
    # PLOT
    fig = plt.figure()
    # Plot ground
    ax = fig.add_subplot(111)
    ax.vlines(pixel_position, 0, match_calipso.calipso.elevation,
              color='k', alpha=1.0)
    # plot cloudsat if we have it
    if match_clsat is not None:
        #  Colors
        colors = []
        for i in range(10):
            colors.append(np.divide([(40 - i) * 6, (40 - i) * 6, (40 - i) * 6], 255.))
        for i in range(10, 30):
            colors.append(np.divide([100 + (i - 10) * 5, 100 + (i - 10) * 5, 100 - (i - 10) * 5], 255.))
        for i in range(30, 100):
            colors.append(np.divide([255, 50, 50], 255.))

        # Plot CloudSat
        base_height = match_clsat.cloudsat.validation_height_base
        top_height = match_clsat.cloudsat.validation_height
        plot_these = top_height > 0
        if RESOLUTION == 5:
            clsat_pos = np.arange(match_clsat.cloudsat.latitude.shape[0])[plot_these]
        else:    
            clsat_pos = pixel_position[match_clsat.cloudsat.calipso_index[plot_these]]
        ax.vlines(clsat_pos,
                  base_height[plot_these], top_height[plot_these], 'm',  # color = colors[nidx], \
                  linestyle='solid', linewidth=1, rasterized=True, label='CPR (CloudSat)')
        title = "%s-CloudSat-CALIOP Cloud Top Heights" % instrument.upper()
    else:
        title = "%s-CALIOP Cloud Top Heights" % instrument.upper()
    #  Plot Caliop
    caliop_label_set = False
    # Plot all 10 calipso layers
    for i in range(10):
        base_ok = caliop_base[:, i]
        top_ok = caliop_height[:, i]
        all_thin = calipso_val_h <= base_ok
        all_thick = np.logical_or(calipso_val_h >= top_ok, calipso_val_h <= 0)
        half_thin = np.logical_and(calipso_val_h > base_ok,
                                   calipso_val_h < top_ok)
        all_thin = np.logical_and(all_thin, top_ok > 0)
        all_thick = np.logical_and(all_thick, top_ok > 0)
        half_thin = np.logical_and(half_thin, top_ok > 0)

        if np.min(top_ok < 0):
            # no more clouds, quit plotting calipso
            break
        if caliop_label_set:
            ax.vlines(pixel_position[all_thick],
                      base_ok[all_thick], top_ok[all_thick], linewidth=0.5,
                      colors="g", linestyle='solid',
                      alpha=0.5, rasterized=True)
        else:
            ax.vlines(pixel_position[all_thick],
                      base_ok[all_thick], top_ok[all_thick],
                      colors="g", linewidth=0.5, linestyle='solid',
                      alpha=0.5, label='caliop', rasterized=True)
            caliop_label_set = True
        ax.vlines(pixel_position[all_thin],
                  base_ok[all_thin], top_ok[all_thin], linewidth=0.5,
                  colors="y", linestyle='solid',
                  alpha=0.5, rasterized=True)
        ax.vlines(pixel_position[half_thin],
                  calipso_val_h[half_thin], top_ok[half_thin], linewidth=0.5,
                  colors="y", linestyle='solid',
                  alpha=0.5, rasterized=True)
        # ax.vlines(pixel_position[half_thin],
        #          base_ok[half_thin], calipso_val_h[half_thin], linewidth=0.5,
        #          colors="g", linestyle='solid',
        #          alpha=0.5, rasterized=True)

    #  Plot Imager
    got_height = imager_ctth_m_above_seasurface >= 0
    ax.plot(pixel_position[got_height], imager_ctth_m_above_seasurface[got_height], 'b+',
            label=instrument.upper(), rasterized=True)
    ax.set_ylim(0, maxheight)
    # plt.show()
    ax.set_title(title)
    ax.set_xlabel("Track Position")
    ax.set_ylabel("Cloud Height (meter)")
    plt.legend(fancybox=True, loc=1, numpoints=4)
    if isinstance(file_type, str):
        filename = "%s/%skm_%s_cloudsat_calipso_%s_clouds_%s.%s" \
            % (plotpath, RESOLUTION, basename, instrument, mode.lower(), file_type)
        fig.savefig(filename, format=file_type)
    else:
        for filetype in file_type:
            filename = "%s/%skm_%s_cloudsat_calipso_%s_clouds_%s.%s" \
                % (plotpath, RESOLUTION, basename, instrument, mode.lower(), filetype)
            fig.savefig(filename, format=filetype) def plot_cal_clsat_imager_satz(match_clsat,
                               match_calipso,
                               plotpath, basename,
                               resolution,
                               file_type='eps',
                               **options):
    """Plot the satellite zenith angle for the track."""
    if 'instrument' in options:
        instrument = options['instrument']
    else:
        instrument = 'imager'
    pixel_position = np.arange(match_calipso.calipso.latitude.shape[0])
    # Plot Satellite Zenith Angle
    fig = plt.figure()
    ax = fig.add_subplot(111)
    if match_clsat is not None:
        use = match_clsat.cloudsat.calipso_index>0
        ax.plot(pixel_position[match_clsat.cloudsat.calipso_index][use],
                match_clsat.imager.satz[use],
                'r+',
                label="%s satz - CloudSat" % instrument.upper())
    ax.set_xlabel("Track Position")
    ax.set_ylabel("satellite zenith angle [deg]")
    ax.set_title("%s SATZ" % instrument.upper())
    ax.plot(match_calipso.imager.satz, "g", label="%s - CALIPSO" % instrument.upper())
    ax.legend(numpoints=4)
    if isinstance(file_type, str):
        fig.savefig("%s/%skm_%s_satz.%s" % (plotpath, resolution, basename, file_type))
    else:
        for filetype in file_type:
            fig.savefig("%s/%skm_%s_satz.%s" % (plotpath, resolution, basename, filetype))
    return def plot_cal_clsat_cwc_imager(match_clsat, elevationcwc, data_okcwc,
                              plotpath, basename, phase, **options):
    if phase == 'IW':
        y1 = match_clsat.cloudsat.RVOD_ice_water_path
        # dataC=match_clsat.cloudsat.RVOD_ice_water_content
        y2 = match_clsat.imager.cpp_iwp
    else:
        y1 = match_clsat.cloudsat.RVOD_liq_water_path
        # dataC=match_clsat.cloudsat.RVOD_liq_water_content
        y2 = match_clsat.imager.cpp_lwp
    pixel_position = np.arange(match_clsat.cloudsat.latitude.shape[0])
    use = np.logical_and(y1 >= 0, y2 >= 0)
    # Findes max value and add 100
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(pixel_position[use], y1[use], "k.", label='RVOD')
    ax.plot(pixel_position[use], y2[use], "b.", label='PPS')
    # ax.set_ylim(0, maxvalue)
    ax.set_xlabel("Track Position")
    ax.set_ylabel("%sP [g/m^2]" % (phase))
    ax.set_title("CloudSat %sP" % (phase))
    filename = "%s/%ikm_%s_cloudsat_rvod_%sP." % (plotpath, RESOLUTION, basename, phase)
    fig.savefig(filename + 'png')
    return def read_truth_imager_match_obj(filename, truth='calipso',
                                read_all=True,
                                read_var=[],
                                skip_var=[]):
    retv = TruthImagerTrackObject(truth=truth)
    h5file = h5py.File(filename, 'r')
    (h5_groups, data_objects) = get_stuff_to_read_from_a_reshaped_file(h5file, retv)
    for group, data_obj in zip(h5_groups, data_objects):
        for dataset in group.keys():
            if dataset in skip_var:
                continue
            if (read_all or dataset in read_var or
                    (len(read_var) == 0 and dataset.data_obj.all_arrays.keys())):
                atrain_match_name = dataset
                if atrain_match_name in ["snow_ice_surface_type"]:
                    atrain_match_name = "nsidc_surface_type"
                setattr(data_obj, atrain_match_name, group[dataset][...])
    retv.diff_sec_1970 = h5file['diff_sec_1970'][...]
    h5file.close()
    return retv class CalipsoObject(DataObject):
    def __init__(self):
        DataObject.__init__(self)
        self.all_arrays = {
            # Normal name = calipso.name.lower()

            # Imager matching needed for all truths:
            'longitude': None,
            'latitude': None,
            'imager_linnum': None,
            'imager_pixnum': None,
            'elevation': None,  # DEM_elevation => elevation in (m)"
            'cloud_fraction': None,
            'validation_height': None,
            'sec_1970': None,
            'minimum_laser_energy_532': None,
            'layer_top_altitude': None,
            'layer_top_temperature': None,
            'layer_top_pressure': None,
            'midlayer_temperature': None,
            'layer_base_altitude': None,
            'layer_base_pressure': None,
            'number_layers_found': None,
            'igbp_surface_type': None,
            'nsidc_surface_type': None,  # V4 renamed from 'snow_ice_surface_type'
            'snow_ice_surface_type': None,
            # 'nsidc_surface_type_texture': None,
            'profile_time_tai': None,  # renamed from "Profile_Time"
            'feature_classification_flags': None,
            'day_night_flag': None,
            'feature_optical_depth_532': None,
            'tropopause_height': None,
            'profile_id': None,

            # If a combination of 5 and 1km data are used for RESOLUTION=1
            # "column_optical_depth_tropospheric_aerosols_1064_5km": None,
            # "column_optical_depth_tropospheric_aerosols_1064": None,
            "column_optical_depth_tropospheric_aerosols_532_5km": None,
            "column_optical_depth_tropospheric_aerosols_532": None,
            "column_optical_depth_aerosols_532_5km": None,
            "column_optical_depth_aerosols_532": None,
            # "column_optical_depth_tropospheric_aerosols_uncertainty_1064_5km": None,
            # "column_optical_depth_tropospheric_aerosols_uncertainty_532_5km": None,
            "column_optical_depth_cloud_532_5km": None,
            # "column_optical_depth_cloud_uncertainty_532_5km": None,
            "feature_optical_depth_532_5km": None,
            "layer_top_altitude_5km": None,
            "layer_top_pressure_5km": None,
            "number_layers_found_5km": None,
            # Variables derived for 5km data
            # Also included if a combination of 5 and 1km data are used for RESOLUTION=1
            'detection_height_5km': None,
            'total_optical_depth_5km': None,
            "feature_optical_depth_532_top_layer_5km": None,
            'cfc_single_shots_1km_from_5km_file': None,
            "average_cloud_top_pressure_single_shots": None,
            "average_cloud_top_pressure_single_shots_5km": None,
            "average_cloud_top_single_shots": None,
            "average_cloud_top_single_shots_5km": None,
            "average_cloud_base_single_shots": None,
            "average_cloud_base_single_shots_5km": None,
            "single_shot_data": None,
            # Variables derived from 5km file to 1kmresolution_
            'cfc_single_shots_1km_from_5km_file': None,

            # From cloudsat:
            'cal_modis_cflag': None,
            'cloudsat_index': None,
        }
import atrain_match.config as config
import os
import numpy as np
import logging
from atrain_match.utils.common import MatchupError
from atrain_match.libs.truth_imager_match import (get_matchups_from_data,
                                                  find_main_cloudproduct_file,
                                                  insert_info_in_filename_or_path,
                                                  add_additional_clousat_calipso_index_vars,
                                                  add_elevation_corrected_imager_ctth)
from atrain_match.truths.cloudsat import (add_validation_ctth_cloudsat,
                                          add_cloudsat_cloud_fraction)
from atrain_match.truths.calipso import (optical_depth_height_filtering,
                                         check_total_optical_depth_and_warn,
                                         add_validation_ctth_calipso,
                                         detection_height_filtering,
                                         set_thin_to_clear_filtering_1km)
from atrain_match.libs.truth_imager_statistics_lib import (calculate_statistics)
from atrain_match.plotting.trajectory_plotting import plot_satellite_trajectory
from atrain_match.plotting.along_track_plotting import (plot_cal_clsat_imager_time_diff,
                                                        plot_cal_clsat_geoprof_imager,
                                                        plot_cal_clsat_imager_satz,
                                                        plot_cal_clsat_cwc_imager)
from atrain_match.matchobject_io import (read_truth_imager_match_obj,
                                         CalipsoObject)
    # TRAJECTORY
    if match_calipso is not None and 1 == 2:
        imlon = match_calipso.imager.longitude.copy()
        imlat = match_calipso.imager.latitude.copy()
        trajectorypath = os.path.join(plotpath, "trajectory_plot")
        if not os.path.exists(trajectorypath):
            os.makedirs(trajectorypath)
        trajectoryname = os.path.join(trajectorypath,
                                      "%skm_%s_trajectory" % (int(config.RESOLUTION),
                                                              values['basename']))
        plot_satellite_trajectory(imlon,
                                  imlat,
                                  trajectoryname,
                                  config.AREA_CONFIG_FILE,
                                  file_type,
                                  **AM_PATHS)

    if (match_calipso is not None):
        # HEIGHT
        plot_cal_clsat_geoprof_imager(match_clsat,
                                      match_calipso,
                                      match_calipso.imager.imager_ctth_m_above_seasurface,
                                      plotpath,
                                      basename,
                                      process_mode,
                                      file_type,
                                      instrument=match_calipso.imager_instrument,
                                      MAXHEIGHT=SETTINGS["MAXHEIGHT"])
        # TIME DIFF SATZ
        plot_cal_clsat_imager_time_diff(match_clsat,
                                        match_calipso,
                                        plotpath, basename,
                                        config.RESOLUTION,
                                        instrument=match_calipso.imager_instrument)
        plot_cal_clsat_imager_satz(match_clsat,
                                   match_calipso,
                                   plotpath, basename,
                                   config.RESOLUTION, file_type,
                                   instrument=match_calipso.imager_instrument)

    if (match_clsat is not None and
            'rvod_liq_water_path' in match_clsat.cloudsat.all_arrays.keys()):

        elevation = np.where(np.less_equal(match_clsat.cloudsat.elevation, 0),
                             -9, match_clsat.cloudsat.elevation)
        data_ok = np.ones(match_clsat.cloudsat.elevation.shape, 'b')

        phase = 'LW'
        plot_cal_clsat_cwc_imager(match_clsat,
                                  elevation,
                                  data_ok,
                                  plotpath, basename,
                                  phase,
                                  instrument=match_clsat.imager_instrument)
        phase = 'IW'
        plot_cal_clsat_cwc_imager(match_clsat,
                                  elevation,
                                  data_ok,
                                  plotpath, basename, phase,
                                  instrument=match_clsat.imager_instrument)


def split_process_mode_and_dnt_part(process_mode_dnt):
    """Slit the process mode to two parts mode and day/night/twilight (DNT)."""
    mode_dnt = process_mode_dnt.split('_')
    if len(mode_dnt) == 1:
        process_mode = process_mode_dnt
        dnt_flag = None
    elif mode_dnt[-1] in ['DAY', 'NIGHT', 'TWILIGHT']:
        process_mode = '_'.join(mode_dnt[0:-1])
        dnt_flag = mode_dnt[-1]
    else:
        process_mode = process_mode_dnt
        dnt_flag = None
    return process_mode, dnt_flag


def process_one_mode(process_mode_dnt, match_calipso, match_clsat, iss_obj, am_obj, sy_obj,
                     min_optical_depth, values, AM_PATHS, SETTINGS, basename):
    """Make plots and statistics for one imager cloudproduct matchup (with one or several truths)."""

    # Get result filename
    process_mode, dnt_flag = split_process_mode_and_dnt_part(process_mode_dnt)
    min_depth_to_file_name = ""
    if process_mode == 'OPTICAL_DEPTH':
        min_depth_to_file_name = "-%.2f" % (min_optical_depth)
    values['mode'] = process_mode_dnt + min_depth_to_file_name
    result_path = insert_info_in_filename_or_path(AM_PATHS['result_dir'],
                                                  values,
                                                  datetime_obj=values['date_time'])
    if not os.path.exists(result_path):
        os.makedirs(result_path)
    result_file = AM_PATHS['result_file'].format(
        resolution=str(config.RESOLUTION),
        basename=values['basename'],
        truth_sat="xxx")
    statfilename = os.path.join(result_path, result_file)
    # Draw plots
    logger.debug("Plotting")
    if process_mode_dnt in SETTINGS['PLOT_MODES']:
        plot_some_figures(match_clsat, match_calipso, values, basename, process_mode,
                          AM_PATHS, SETTINGS, am_obj=am_obj)
    # Calculate Statistics
    logger.debug("Calculating statistics")
    calculate_statistics(process_mode, statfilename, match_calipso, match_clsat,
                         iss_obj, am_obj, sy_obj, SETTINGS, dnt_flag)


def run(cross, run_modes, AM_PATHS, SETTINGS, reprocess=False):
    """The main work horse."""
    logger.info("Case: %s", str(cross))
    # sensor = INSTRUMENT.get(cross.satellite1.lower(), 'imager')

    if (not SETTINGS['USE_CMA_FOR_CFC_STATISTICS'] and
        not SETTINGS['USE_CT_FOR_CFC_STATISTICS'] and
            not SETTINGS['USE_CMAPROB_FOR_CFC_STATISTICS']):
        logger.error(
            "\n  --------------------------------- "
            "\n\tSet one of USE_*_FOR_CFC_STATISTICS=True in config.py!"
            "\n  --------------------------------- ")

        raise MatchupError("Configure problems, see messages above.")