def rectwv_coeff_from_arc_image(reduced_image,
                                bound_param,
                                lines_catalog,
                                args_nbrightlines=None,
                                args_ymargin_bb=2,
                                args_remove_sp_background=True,
                                args_times_sigma_threshold=10,
                                args_order_fmap=2,
                                args_sigma_gaussian_filtering=2,
                                args_margin_npix=50,
                                args_poldeg_initial=3,
                                args_poldeg_refined=5,
                                args_interactive=False,
                                args_threshold_wv=0,
                                args_ylogscale=False,
                                args_pdf=None,
                                args_geometry=(0,0,640,480),
                                debugplot=0):
    """Evaluate rect.+wavecal. coefficients from arc image

    Parameters
    ----------
    reduced_image : HDUList object
        Image with preliminary basic reduction: bpm, bias, dark and
        flatfield.
    bound_param : RefinedBoundaryModelParam instance
        Refined boundary model.
    lines_catalog : Numpy array
        2D numpy array with the contents of the master file with the
        expected arc line wavelengths.
    args_nbrightlines : int
        TBD
    args_ymargin_bb : int
        TBD
    args_remove_sp_background : bool
        TBD
    args_times_sigma_threshold : float
        TBD
    args_order_fmap : int
        TBD
    args_sigma_gaussian_filtering : float
        TBD
    args_margin_npix : int
        TBD
    args_poldeg_initial : int
        TBD
    args_poldeg_refined : int
        TBD
    args_interactive : bool
        TBD
    args_threshold_wv : float
        TBD
    args_ylogscale : bool
        TBD
    args_pdf : TBD
    args_geometry : TBD
    debugplot : int
            Debugging level for messages and plots. For details see
            'numina.array.display.pause_debugplot.py'.

    Returns
    -------
    rectwv_coeff : RectWaveCoeff instance
        Rectification and wavelength calibration coefficients for the
        particular CSU configuration of the input arc image.
    reduced_55sp : HDUList object
        Image with 55 spectra corresponding to the median spectrum for
        each slitlet, employed to derived the wavelength calibration
        polynomial.

    """

    logger = logging.getLogger(__name__)

    # protections
    if args_interactive and args_pdf is not None:
        logger.error('--interactive and --pdf are incompatible options')
        raise ValueError('--interactive and --pdf are incompatible options')

    # header and data array
    header = reduced_image[0].header
    image2d = reduced_image[0].data

    # check grism and filter
    filter_name = header['filter']
    logger.info('Filter: ' + filter_name)
    if filter_name != bound_param.tags['filter']:
        raise ValueError('Filter name does not match!')
    grism_name = header['grism']
    logger.info('Grism: ' + grism_name)
    if grism_name != bound_param.tags['grism']:
        raise ValueError('Grism name does not match!')

    # read the CSU configuration from the image header
    csu_conf = CsuConfiguration.define_from_header(header)
    logger.debug(csu_conf)

    # read the DTU configuration from the image header
    dtu_conf = DtuConfiguration.define_from_header(header)
    logger.debug(dtu_conf)

    # set boundary parameters
    parmodel = bound_param.meta_info['parmodel']
    params = bound_params_from_dict(bound_param.__getstate__())
    if abs(debugplot) >= 10:
        print('-' * 83)
        print('* FITTED BOUND PARAMETERS')
        params.pretty_print()
        pause_debugplot(debugplot)

    # determine parameters according to grism+filter combination
    wv_parameters = set_wv_parameters(filter_name, grism_name)
    islitlet_min = wv_parameters['islitlet_min']
    islitlet_max = wv_parameters['islitlet_max']
    if args_nbrightlines is None:
        nbrightlines = wv_parameters['nbrightlines']
    else:
        nbrightlines = [int(idum) for idum in args_nbrightlines.split(',')]
    poly_crval1_linear = wv_parameters['poly_crval1_linear']
    poly_cdelt1_linear = wv_parameters['poly_cdelt1_linear']
    wvmin_expected = wv_parameters['wvmin_expected']
    wvmax_expected = wv_parameters['wvmax_expected']
    wvmin_useful = wv_parameters['wvmin_useful']
    wvmax_useful = wv_parameters['wvmax_useful']

    # list of slitlets to be computed
    logger.info('list_slitlets: [' + str(islitlet_min) + ',... ' +
                str(islitlet_max) + ']')

    # read master arc line wavelengths (only brightest lines)
    wv_master = read_wv_master_from_array(
        master_table=lines_catalog, lines='brightest', debugplot=debugplot
    )

    # read master arc line wavelengths (whole data set)
    wv_master_all = read_wv_master_from_array(
        master_table=lines_catalog, lines='all', debugplot=debugplot
    )

    # check that the arc lines in the master file are properly sorted
    # in ascending order
    for i in range(len(wv_master_all) - 1):
        if wv_master_all[i] >= wv_master_all[i + 1]:
            logger.error('>>> wavelengths: ' +
                         str(wv_master_all[i]) + '  ' +
                         str(wv_master_all[i+1]))
            raise ValueError('Arc lines are not sorted in master file')

    # ---

    image2d_55sp = np.zeros((EMIR_NBARS, EMIR_NAXIS1))

    # compute rectification transformation and wavelength calibration
    # polynomials

    measured_slitlets = []

    cout = '0'
    for islitlet in range(1, EMIR_NBARS + 1):

        if islitlet_min <= islitlet <= islitlet_max:

            # define Slitlet2dArc object
            slt = Slitlet2dArc(
                islitlet=islitlet,
                csu_conf=csu_conf,
                ymargin_bb=args_ymargin_bb,
                params=params,
                parmodel=parmodel,
                debugplot=debugplot
            )

            # extract 2D image corresponding to the selected slitlet, clipping
            # the image beyond the unrectified slitlet (in order to isolate
            # the arc lines of the current slitlet; otherwise there are
            # problems with arc lines from neighbour slitlets)
            image2d_tmp = select_unrectified_slitlet(
                image2d=image2d,
                islitlet=islitlet,
                csu_bar_slit_center=csu_conf.csu_bar_slit_center(islitlet),
                params=params,
                parmodel=parmodel,
                maskonly=False
            )
            slitlet2d = slt.extract_slitlet2d(image2d_tmp)

            # subtract smooth background computed as follows:
            # - median collapsed spectrum of the whole slitlet2d
            # - independent median filtering of the previous spectrum in the
            #   two halves in the spectral direction
            if args_remove_sp_background:
                spmedian = np.median(slitlet2d, axis=0)
                naxis1_tmp = spmedian.shape[0]
                jmidpoint = naxis1_tmp // 2
                sp1 = medfilt(spmedian[:jmidpoint], [201])
                sp2 = medfilt(spmedian[jmidpoint:], [201])
                spbackground = np.concatenate((sp1, sp2))
                slitlet2d -= spbackground

            # locate unknown arc lines
            slt.locate_unknown_arc_lines(
                slitlet2d=slitlet2d,
                times_sigma_threshold=args_times_sigma_threshold)

            # continue working with current slitlet only if arc lines have
            # been detected
            if slt.list_arc_lines is not None:

                # compute intersections between spectrum trails and arc lines
                slt.xy_spectrail_arc_intersections(slitlet2d=slitlet2d)

                # compute rectification transformation
                slt.estimate_tt_to_rectify(order=args_order_fmap,
                                           slitlet2d=slitlet2d)

                # rectify image
                slitlet2d_rect = slt.rectify(slitlet2d,
                                             resampling=2,
                                             transformation=1)

                # median spectrum and line peaks from rectified image
                sp_median, fxpeaks = slt.median_spectrum_from_rectified_image(
                    slitlet2d_rect,
                    sigma_gaussian_filtering=args_sigma_gaussian_filtering,
                    nwinwidth_initial=5,
                    nwinwidth_refined=5,
                    times_sigma_threshold=5,
                    npix_avoid_border=6,
                    nbrightlines=nbrightlines
                )

                image2d_55sp[islitlet - 1, :] = sp_median

                # determine expected wavelength limits prior to the wavelength
                # calibration
                csu_bar_slit_center = csu_conf.csu_bar_slit_center(islitlet)
                crval1_linear = poly_crval1_linear(csu_bar_slit_center)
                cdelt1_linear = poly_cdelt1_linear(csu_bar_slit_center)
                expected_wvmin = crval1_linear - \
                                 args_margin_npix * cdelt1_linear
                naxis1_linear = sp_median.shape[0]
                crvaln_linear = crval1_linear + \
                                (naxis1_linear - 1) * cdelt1_linear
                expected_wvmax = crvaln_linear + \
                                 args_margin_npix * cdelt1_linear
                # override previous estimates when necessary
                if wvmin_expected is not None:
                    expected_wvmin = wvmin_expected
                if wvmax_expected is not None:
                    expected_wvmax = wvmax_expected

                # clip initial master arc line list with bright lines to
                # the expected wavelength range
                lok1 = expected_wvmin <= wv_master
                lok2 = wv_master <= expected_wvmax
                lok = lok1 * lok2
                wv_master_eff = wv_master[lok]

                # perform initial wavelength calibration
                solution_wv = wvcal_spectrum(
                    sp=sp_median,
                    fxpeaks=fxpeaks,
                    poly_degree_wfit=args_poldeg_initial,
                    wv_master=wv_master_eff,
                    wv_ini_search=expected_wvmin,
                    wv_end_search=expected_wvmax,
                    wvmin_useful=wvmin_useful,
                    wvmax_useful=wvmax_useful,
                    geometry=args_geometry,
                    debugplot=slt.debugplot
                )
                # store initial wavelength calibration polynomial in current
                # slitlet instance
                slt.wpoly = np.polynomial.Polynomial(solution_wv.coeff)
                pause_debugplot(debugplot)

                # clip initial master arc line list with all the lines to
                # the expected wavelength range
                lok1 = expected_wvmin <= wv_master_all
                lok2 = wv_master_all <= expected_wvmax
                lok = lok1 * lok2
                wv_master_all_eff = wv_master_all[lok]

                # clip master arc line list to useful region
                if wvmin_useful is not None:
                    lok = wvmin_useful <= wv_master_all_eff
                    wv_master_all_eff  = wv_master_all_eff[lok]
                if wvmax_useful is not None:
                    lok = wv_master_all_eff <= wvmax_useful
                    wv_master_all_eff  = wv_master_all_eff[lok]

                # refine wavelength calibration
                if args_poldeg_refined > 0:
                    plottitle = '[slitlet#{}, refined]'.format(islitlet)
                    poly_refined, yres_summary = refine_arccalibration(
                        sp=sp_median,
                        poly_initial=slt.wpoly,
                        wv_master=wv_master_all_eff,
                        poldeg=args_poldeg_refined,
                        ntimes_match_wv=1,
                        interactive=args_interactive,
                        threshold=args_threshold_wv,
                        plottitle=plottitle,
                        ylogscale=args_ylogscale,
                        geometry=args_geometry,
                        pdf=args_pdf,
                        debugplot=slt.debugplot
                    )
                    # store refined wavelength calibration polynomial in
                    # current slitlet instance
                    slt.wpoly = poly_refined

                # compute approximate linear values for CRVAL1 and CDELT1
                naxis1_linear = sp_median.shape[0]
                crmin1_linear = slt.wpoly(1)
                crmax1_linear = slt.wpoly(naxis1_linear)
                slt.crval1_linear = crmin1_linear
                slt.cdelt1_linear = \
                    (crmax1_linear - crmin1_linear) / (naxis1_linear - 1)

                # check that the trimming of wv_master and wv_master_all has
                # preserved the wavelength range [crmin1_linear, crmax1_linear]
                if crmin1_linear < expected_wvmin:
                    logger.warning(">>> islitlet: " +str(islitlet))
                    logger.warning("expected_wvmin: " + str(expected_wvmin))
                    logger.warning("crmin1_linear.: " + str(crmin1_linear))
                    logger.warning("WARNING: Unexpected crmin1_linear < "
                                   "expected_wvmin")
                if crmax1_linear > expected_wvmax:
                    logger.warning(">>> islitlet: " +str(islitlet))
                    logger.warning("expected_wvmax: " + str(expected_wvmax))
                    logger.warning("crmax1_linear.: " + str(crmax1_linear))
                    logger.warning("WARNING: Unexpected crmax1_linear > "
                                   "expected_wvmax")

                cout += '.'

            else:

                cout += 'x'

            if islitlet % 10 == 0:
                if cout != 'x':
                    cout = str(islitlet // 10)

            if debugplot != 0:
                pause_debugplot(debugplot)

        else:

            # define Slitlet2dArc object
            slt = Slitlet2dArc(
                islitlet=islitlet,
                csu_conf=csu_conf,
                ymargin_bb=args_ymargin_bb,
                params=None,
                parmodel=None,
                debugplot=debugplot
            )

            cout += 'i'

        # store current slitlet in list of measured slitlets
        measured_slitlets.append(slt)

        logger.info(cout)

    # ---

    # generate FITS file structure with 55 spectra corresponding to the
    # median spectrum for each slitlet
    reduced_55sp = fits.PrimaryHDU(data=image2d_55sp)
    reduced_55sp.header['crpix1'] = (0.0, 'reference pixel')
    reduced_55sp.header['crval1'] = (0.0, 'central value at crpix2')
    reduced_55sp.header['cdelt1'] = (1.0, 'increment')
    reduced_55sp.header['ctype1'] = 'PIXEL'
    reduced_55sp.header['cunit1'] = ('Pixel', 'units along axis2')
    reduced_55sp.header['crpix2'] = (0.0, 'reference pixel')
    reduced_55sp.header['crval2'] = (0.0, 'central value at crpix2')
    reduced_55sp.header['cdelt2'] = (1.0, 'increment')
    reduced_55sp.header['ctype2'] = 'PIXEL'
    reduced_55sp.header['cunit2'] = ('Pixel', 'units along axis2')

    # ---

    # Generate structure to store intermediate results
    outdict = {}
    outdict['instrument'] = 'EMIR'
    outdict['meta_info'] = {}
    outdict['meta_info']['creation_date'] = datetime.now().isoformat()
    outdict['meta_info']['description'] = \
        'computation of rectification and wavelength calibration polynomial ' \
        'coefficients for a particular CSU configuration'
    outdict['meta_info']['recipe_name'] = 'undefined'
    outdict['meta_info']['origin'] = {}
    outdict['meta_info']['origin']['bound_param_uuid'] = \
        bound_param.uuid
    outdict['meta_info']['origin']['arc_image_uuid'] = 'undefined'
    outdict['tags'] = {}
    outdict['tags']['grism'] = grism_name
    outdict['tags']['filter'] = filter_name
    outdict['tags']['islitlet_min'] = islitlet_min
    outdict['tags']['islitlet_max'] = islitlet_max
    outdict['dtu_configuration'] = dtu_conf.outdict()
    outdict['uuid'] = str(uuid4())
    outdict['contents'] = {}

    missing_slitlets = []
    for slt in measured_slitlets:

        islitlet = slt.islitlet

        if islitlet_min <= islitlet <= islitlet_max:

            # avoid error when creating a python list of coefficients from
            # numpy polynomials when the polynomials do not exist (note that
            # the JSON format doesn't handle numpy arrays and such arrays must
            # be transformed into native python lists)
            if slt.wpoly is None:
                wpoly_coeff = None
            else:
                wpoly_coeff = slt.wpoly.coef.tolist()
            if slt.wpoly_longslit_model is None:
                wpoly_coeff_longslit_model = None
            else:
                wpoly_coeff_longslit_model = \
                    slt.wpoly_longslit_model.coef.tolist()

            # avoid similar error when creating a python list of coefficients
            # when the numpy array does not exist; note that this problem
            # does not happen with tt?_aij_longslit_model and
            # tt?_bij_longslit_model because the latter have already been
            # created as native python lists
            if slt.ttd_aij is None:
                ttd_aij = None
            else:
                ttd_aij = slt.ttd_aij.tolist()
            if slt.ttd_bij is None:
                ttd_bij = None
            else:
                ttd_bij = slt.ttd_bij.tolist()
            if slt.tti_aij is None:
                tti_aij = None
            else:
                tti_aij = slt.tti_aij.tolist()
            if slt.tti_bij is None:
                tti_bij = None
            else:
                tti_bij = slt.tti_bij.tolist()

            # creating temporary dictionary with the information corresponding
            # to the current slitlett that will be saved in the JSON file
            tmp_dict = {
                'csu_bar_left': slt.csu_bar_left,
                'csu_bar_right': slt.csu_bar_right,
                'csu_bar_slit_center': slt.csu_bar_slit_center,
                'csu_bar_slit_width': slt.csu_bar_slit_width,
                'x0_reference': slt.x0_reference,
                'y0_reference_lower': slt.y0_reference_lower,
                'y0_reference_middle': slt.y0_reference_middle,
                'y0_reference_upper': slt.y0_reference_upper,
                'y0_reference_lower_expected':
                    slt.y0_reference_lower_expected,
                'y0_reference_middle_expected':
                    slt.y0_reference_middle_expected,
                'y0_reference_upper_expected':
                    slt.y0_reference_upper_expected,
                'y0_frontier_lower': slt.y0_frontier_lower,
                'y0_frontier_upper': slt.y0_frontier_upper,
                'y0_frontier_lower_expected': slt.y0_frontier_lower_expected,
                'y0_frontier_upper_expected': slt.y0_frontier_upper_expected,
                'corr_yrect_a': slt.corr_yrect_a,
                'corr_yrect_b': slt.corr_yrect_b,
                'min_row_rectified': slt.min_row_rectified,
                'max_row_rectified': slt.max_row_rectified,
                'ymargin_bb': slt.ymargin_bb,
                'bb_nc1_orig': slt.bb_nc1_orig,
                'bb_nc2_orig': slt.bb_nc2_orig,
                'bb_ns1_orig': slt.bb_ns1_orig,
                'bb_ns2_orig': slt.bb_ns2_orig,
                'spectrail': {
                    'poly_coef_lower':
                        slt.list_spectrails[
                            slt.i_lower_spectrail].poly_funct.coef.tolist(),
                    'poly_coef_middle':
                        slt.list_spectrails[
                            slt.i_middle_spectrail].poly_funct.coef.tolist(),
                    'poly_coef_upper':
                        slt.list_spectrails[
                            slt.i_upper_spectrail].poly_funct.coef.tolist(),
                },
                'frontier': {
                    'poly_coef_lower':
                        slt.list_frontiers[0].poly_funct.coef.tolist(),
                    'poly_coef_upper':
                        slt.list_frontiers[1].poly_funct.coef.tolist(),
                },
                'ttd_order': slt.ttd_order,
                'ttd_aij': ttd_aij,
                'ttd_bij': ttd_bij,
                'tti_aij': tti_aij,
                'tti_bij': tti_bij,
                'ttd_order_longslit_model': slt.ttd_order_longslit_model,
                'ttd_aij_longslit_model': slt.ttd_aij_longslit_model,
                'ttd_bij_longslit_model': slt.ttd_bij_longslit_model,
                'tti_aij_longslit_model': slt.tti_aij_longslit_model,
                'tti_bij_longslit_model': slt.tti_bij_longslit_model,
                'wpoly_coeff': wpoly_coeff,
                'wpoly_coeff_longslit_model': wpoly_coeff_longslit_model,
                'crval1_linear': slt.crval1_linear,
                'cdelt1_linear': slt.cdelt1_linear
            }
        else:
            missing_slitlets.append(islitlet)
            tmp_dict = {
                'csu_bar_left': slt.csu_bar_left,
                'csu_bar_right': slt.csu_bar_right,
                'csu_bar_slit_center': slt.csu_bar_slit_center,
                'csu_bar_slit_width': slt.csu_bar_slit_width,
                'x0_reference': slt.x0_reference,
                'y0_frontier_lower_expected': slt.y0_frontier_lower_expected,
                'y0_frontier_upper_expected': slt.y0_frontier_upper_expected
            }
        slitlet_label = "slitlet" + str(islitlet).zfill(2)
        outdict['contents'][slitlet_label] = tmp_dict

    # ---

    # OBSOLETE
    '''
    # save JSON file needed to compute the MOS model
    with open(args.out_json.name, 'w') as fstream:
        json.dump(outdict, fstream, indent=2, sort_keys=True)
        print('>>> Saving file ' + args.out_json.name)
    '''

    # ---

    # Create object of type RectWaveCoeff with coefficients for
    # rectification and wavelength calibration
    rectwv_coeff = RectWaveCoeff(instrument='EMIR')
    rectwv_coeff.quality_control = numina.types.qc.QC.GOOD
    rectwv_coeff.tags['grism'] = grism_name
    rectwv_coeff.tags['filter'] = filter_name
    rectwv_coeff.meta_info['origin']['bound_param'] = \
        'uuid' + bound_param.uuid
    rectwv_coeff.meta_info['dtu_configuration'] = outdict['dtu_configuration']
    rectwv_coeff.total_slitlets = EMIR_NBARS
    rectwv_coeff.missing_slitlets = missing_slitlets
    for i in range(EMIR_NBARS):
        islitlet = i + 1
        dumdict = {'islitlet': islitlet}
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        if cslitlet in outdict['contents']:
            dumdict.update(outdict['contents'][cslitlet])
        else:
            raise ValueError("Unexpected error")
        rectwv_coeff.contents.append(dumdict)
    # debugging __getstate__ and __setstate__
    # rectwv_coeff.writeto(args.out_json.name)
    # print('>>> Saving file ' + args.out_json.name)
    # check_setstate_getstate(rectwv_coeff, args.out_json.name)
    logger.info('Generating RectWaveCoeff object with uuid=' +
                rectwv_coeff.uuid)

    return rectwv_coeff, reduced_55sp
Evaluate rect.+wavecal. coefficients from arc image

    Parameters
    ----------
    reduced_image : HDUList object
        Image with preliminary basic reduction: bpm, bias, dark and
        flatfield.
    bound_param : RefinedBoundaryModelParam instance
        Refined boundary model.
    lines_catalog : Numpy array
        2D numpy array with the contents of the master file with the
        expected arc line wavelengths.
    args_nbrightlines : int
        TBD
    args_ymargin_bb : int
        TBD
    args_remove_sp_background : bool
        TBD
    args_times_sigma_threshold : float
        TBD
    args_order_fmap : int
        TBD
    args_sigma_gaussian_filtering : float
        TBD
    args_margin_npix : int
        TBD
    args_poldeg_initial : int
        TBD
    args_poldeg_refined : int
        TBD
    args_interactive : bool
        TBD
    args_threshold_wv : float
        TBD
    args_ylogscale : bool
        TBD
    args_pdf : TBD
    args_geometry : TBD
    debugplot : int
            Debugging level for messages and plots. For details see
            'numina.array.display.pause_debugplot.py'.

    Returns
    -------
    rectwv_coeff : RectWaveCoeff instance
        Rectification and wavelength calibration coefficients for the
        particular CSU configuration of the input arc image.
    reduced_55sp : HDUList object
        Image with 55 spectra corresponding to the median spectrum for
        each slitlet, employed to derived the wavelength calibration
        polynomial.