def plot_peak_altitude(self, ax=None, true_color='k',
            apparent_color='grey'):
        if ax is None:
            ax = plt.gca()
        plt.sca(ax)

        for d in self.digitization_list:
            if d.is_invertible():
                if d.altitude.size == 0:
                    try:
                        d.invert(substitute_fp=ais_code.ne_to_fp(4.))
                    except BaseException as e:
                        print(e)
                        continue

                if apparent_color:
                    plt.plot(d.time, d.altitude[-1], marker='.',
                        ms=self.marker_size, color=apparent_color)
                alt = mex.iau_pgr_alt_lat_lon_position(float(d.time))[0]
                plt.plot(d.time,
                alt - d.traced_delay[-1] * ais_code.speed_of_light_kms / 2.,
                        marker='.', color=true_color,
                        ms=self.marker_size)
        celsius.ylabel(r'$h_{max} / km$')
        plt.ylim(0o1, 249)

 def plot_r(self, ax=None, label=True, fmt='k-', **kwargs):
     if ax is None:
         ax = plt.gca()
     plt.sca(ax)
     self.generate_position()
     plt.plot(self.t, self.iau_pos[0] / mex.mars_mean_radius_km, fmt, **kwargs)
     celsius.ylabel(r'$r / R_M$')

    def plot_timeseries(self, ax=None, vmin=None, vmax=None,
            colorbar=False, label=True):

        if vmin is None:
            vmin = self.vmin

        if vmax is None:
            vmax = self.vmax

        if ax is None:
            ax = plt.gca()
        plt.sca(ax)
        plt.cla()
        plt.imshow(self.tser_arr[::-1,:], vmin=vmin, vmax=vmax,
            interpolation='Nearest', extent=self.extent, origin='upper',aspect='auto')
        plt.xlim(self.extent[0], self.extent[1])
        plt.ylim(self.extent[2], self.extent[3])
        # plt.vlines(self.ionogram_list[0].time, self.extent[2], self.extent[3], 'r')
        if label:
           celsius.ylabel('f / MHz')

        if colorbar:
            old_ax = plt.gca()
            plt.colorbar(
                    cax = celsius.make_colorbar_cax(), ticks=self.cbar_ticks
                ).set_label(r"$Log_{10} V^2 m^{-2} Hz^{-1}$")
            plt.sca(old_ax)

def plot_matches_group(g,how='chi',ntop=8):
    """
    Plot matches from h5 group

    Pulls out the relavent arrays from a group and runs plot_matches
    
    Parameters
    ----------
    g : h5 group containing the following datasets
        - arr : DataSet with spectra
        - smres : DataSet with specmatch results
    
    """
    smres = pd.DataFrame(g['smres'][:])
    smres = results.smres_add_chi(smres)

    lspec = g['lspec'][:]
    smres.index = np.arange(len(smres))

    if how=='chi':
        smresbest = smres.sort_values(by='chi')
    elif how=='close':
        targname = smio.kbc_query(smres.targobs[0])['name']
        tpar = dict(lib.ix[targname])
        smres['close'] = close(tpar,smres)
        smresbest = smres.sort_values(by='close')
    smresbest = smresbest.iloc[:ntop]

    plot_matches(smresbest,lspec)
    plt.sca(plt.gcf().get_axes()[0])

def orbit_plots(orbit_list=[8020, 8021, 8022, 8023, 8024, 8025], resolution=200, ax=None):

    if ax is None:
        ax = plt.gca()
    plt.sca(ax)

    orbits = {}
    props = dict(marker="None", hold=True, mec=None, linestyle="-", markeredgewidth=0.0)
    # props = dict(marker='o', hold=True,mec='None',line style='None',markeredgewidth=0.0)

    for orbit in orbit_list:
        orbit_t = mex.orbits[orbit]
        print(orbit, orbit_t.start, orbit_t.finish)
        t = np.linspace(orbit_t.start, orbit_t.finish, resolution)
        pos = mex.iau_pgr_alt_lat_lon_position(t)
        # orbits[orbit] = dict(lat = np.rad2deg(pos[2]), lon = np.rad2deg(np.unwrap(pos[1])),
        #                         alt=pos[0] - mex.mars_mean_radius_km, t=t)
        orbits[orbit] = dict(lat=pos[1], lon=np.rad2deg(np.unwrap(np.deg2rad(pos[2]))), alt=pos[0], t=t)

        ll = plt.plot(orbits[orbit]["lon"], orbits[orbit]["lat"], label=str(orbit), **props)

        plt.plot(
            orbits[orbit]["lon"] + 360.0, orbits[orbit]["lat"], label="_nolegend_", color=ll[0].get_color(), **props
        )
        plt.plot(
            orbits[orbit]["lon"] - 360.0, orbits[orbit]["lat"], label="_nolegend_", color=ll[0].get_color(), **props
        )

    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.xlabel("Longitude / deg")
    plt.ylabel("Latitude / deg")
    plt.legend(loc="center left", bbox_to_anchor=(1.01, 0.5), numpoints=1, title="Orbit")
    ax.xaxis.set_major_locator(celsius.CircularLocator())
    ax.yaxis.set_major_locator(celsius.CircularLocator())

    def plot_acf_ccf(axacf,axccf,wlo,whi):
        b = (w > wlo) & (w < whi)
        plt.sca(axccf)
        dvmax = plot_ccf(w[b], tspec[b], mspec[b])

        plt.sca(axacf)
        thwhm,mhwhm = plot_acf(w[b], tspec[b], mspec[b])
        return dvmax,thwhm,mhwhm

 def setlims(d):
     plt.sca(d['ax'])
     xk = d['xk']
     yk = d['yk']
     if limd.keys().count(xk)==1:
         plt.xlim(*limd[xk])
     if limd.keys().count(yk)==1:
         plt.ylim(*limd[yk])

 def plot_altitude(self, ax=None, label=True, fmt='k-', **kwargs):
     if ax is None:
         ax = plt.gca()
     plt.sca(ax)
     self.generate_position()
     plt.plot(self.t, self.iau_pos[0] - mex.mars_mean_radius_km, fmt, **kwargs)
     if label:
         celsius.ylabel('h / km')

 def plot_sza(self, ax=None, label=True, fmt='k-', **kwargs):
     if ax is None:
         ax = plt.gca()
     plt.sca(ax)
     self.make_axis_circular(ax)
     self.generate_position()
     plt.plot(self.t, self.sza, fmt, **kwargs)
     if label:
         celsius.ylabel(r'$SZA / deg$')
     plt.ylim(0., 180.)

 def flipax(d):
     plt.sca(d['ax'])
     xk = d['xk']
     yk = d['yk']
     if flipd.keys().count(xk)==1:
         if flipd[xk]:
             flip('x')
     if flipd.keys().count(yk)==1:
         if flipd[yk]:
             flip('y')

 def plot_lat(self, ax=None, label=True, fmt='k-', **kwargs):
     if ax is None:
         ax = plt.gca()
     plt.sca(ax)
     self.make_axis_circular(ax)
     self.generate_position()
     plt.plot(self.t, self.iau_pos[1], fmt, **kwargs)
     if label:
         celsius.ylabel(r'$\lambda$')
     plt.ylim(-90., 90.)

    def plot_mod_b(self, fmt='k.', ax=None,
                    field_model=True, errors=True, field_color='blue',
                    br=True, t_offset=0., label=True, **kwargs):
        if ax is None:
            ax = plt.gca()
        plt.sca(ax)

        sub = [d for d in self.digitization_list if np.isfinite(d.td_cyclotron)]
        if len(sub) == 0:
            print("No digitizations with marked cyclotron frequency lines")
            return

        t = np.array([d.time for d in sub])
        b = np.array([d.td_cyclotron for d in sub])
        e = np.array([d.td_cyclotron_error for d in sub])
        # print b
        # print e

        b, e = ais_code.td_to_modb(b, e)
        b *= 1.E9
        e *= 1.E9

        if errors:
            for tt,bb,ee in zip(t,b,e):
                plt.plot((tt,tt),(bb+ee,bb-ee),
                                color='lightgrey',linestyle='solid',marker='None')
                plt.plot(tt,bb,fmt,ms=self.marker_size, **kwargs)
            # plt.errorbar(t, b, e, fmt=fmt, ms=self.marker_size, **kwargs)
        else:
            plt.plot(t, b, fmt, ms=self.marker_size, **kwargs)

        if field_model:
            self.generate_position()

            if field_color is None: field_color = fmt[0]
            # b = self.quick_field_model(self.t)
            self._computed_field_model = self.field_model(self.iau_pos)
            bmag = np.sqrt(np.sum(self._computed_field_model**2., 0))
            plt.plot(self.t - t_offset,
                        bmag,
                        color=field_color, ls='-')
            if br:
                plt.plot(self.t - t_offset,
                    self._computed_field_model[0], 'r-')
                plt.plot(self.t - t_offset,
                    -1. * self._computed_field_model[0], 'r', ls='dashed')

            model_at_value = np.interp(t, self.t, bmag)
            inx = (model_at_value > 100.) & ((b / model_at_value) < 0.75)
            plt.plot(t[inx], b[inx], 'ro', mec='r', mfc='none', ms=5., mew=1.2)


        if label:
            celsius.ylabel(r'$\mathrm{|B|/nT}$')
        plt.ylim(0., 200)

def test_axes_stuff():
    fig = DJAPage(ratios=[2.,1., 0.63])
    plt.sca(fig.top_axes)
    x = np.arange(360.)
    y = np.sin(x * np.pi/180. + 2.)
    plt.plot(x, y, 'r-')
    plt.sca(fig.bottom_axes)
    mb = add_labelled_bar(x,x)
    fig.register_new_axis(mb)
    plt.xlim(0., 360.)
    plt.show()

 def plot_aspera_els(self, ax=None, **kwargs):
     if self.verbose: print('PLOT_ASPERA_ELS:')
     if ax is None:
         ax = plt.gca()
     else:
         plt.sca(ax)
     try:
         mex.aspera.plot_els_spectra(self.extent[0], self.extent[1],
             ax=ax, verbose=self.verbose, **kwargs)
     except Exception as e:
         print(e)

def velocityshift(wav, flux, ref_wav, ref_flux, plot=False):
    """
    Find the velocity shift between two spectra. 

    Args:
        wav (array): Wavelength array.
        flux (array): Continuum-normalized spectrum. 
        ref_wav (array): 
        ref_flux (array): 
    
    Returns:
        vmax (float): Velocity of the cross-correlation peak. Positive velocity
            means that observed spectrum is red-shifted with respect to the 
            reference spectrum.
         corrmax (float): Peak cross-correlation amplitude.
    """
    nwav = flux.size

    # Build spline object for resampling the model spectrum
    ref_spline = InterpolatedUnivariateSpline(ref_wav, ref_flux)

    # Convert target and model spectra to constant log-lambda scale
    wav, flux, dvel = loglambda(wav, flux)
    ref_flux = ref_spline(wav)

    # Perform cross-correlation, and use quadratic interpolation to
    # find the velocity value that maximizes the cross-correlation
    # amplitude. If `lag` is negative, the observed spectrum need to
    # be blue-shifted in order to line up with the observed
    # spectrum. Thus, to put the spectra on the same scale, the
    # observed spectrum must be red-shifted, i.e. vmax is positive.
    flux-=np.mean(flux)
    ref_flux-=np.mean(ref_flux)
    lag = np.arange(-nwav + 1, nwav) 
    dvel = -1.0 * lag * dvel
    corr = np.correlate(ref_flux, flux, mode='full')
    vmax, corrmax = quadratic_max(dvel, corr)

    if plot:
        from matplotlib import pylab as plt
        fig,axL = plt.subplots(ncols=2)
        plt.sca(axL[0])
        plt.plot(wav,ref_flux)
        plt.plot(wav,flux)
        plt.sca(axL[1])
        vrange = (-100,100)
        b = (dvel > vrange[0]) & (dvel < vrange[1])
        plt.plot(dvel[b],corr[b])
        plt.plot([vmax],[corrmax],'o',label='Cross-correlation Peak')
        fig.set_tight_layout(True)
        plt.draw()
        plt.show()

    return vmax, corrmax

 def wrapped_f(*args):
     for ax in axL:                
         plt.sca(ax)
         f(*args)
     xl = plt.xlim()
     start = xl[0]
     step = (xl[1]-xl[0]) / float(nax)
     for ax in axL:                
         plt.sca(ax)
         plt.xlim(start,start+step)
         start+=step

 def plotline(x):
     ax = x['ax']
     xk = x['xk']
     plt.sca(ax)
     trans = blended_transform_factory(ax.transData,ax.transAxes)
     plt.axvline(smpar[xk],ls='--')
     plt.text(smpar[xk],0.9,'SM',transform=trans)
     
     if libpar is not None:
         plt.axvline(libpar[xk])
         plt.text(libpar[xk],0.9,'LIB',transform=trans)

def plot_frame(dfplot,df,**kwargs):
    """
    Plot DataFrame
    
    Helper function for panels.
    """
    for i in dfplot.index:
        x = dfplot.ix[i]
        plt.sca(x['ax'])
        plt.plot(df[x['xk']],df[x['yk']],**kwargs)
        plt.xlabel(x['xL'])
        plt.ylabel(x['yL'])

    def plot_frequency_altitude(self, f=2.0, ax=None, median_filter=False,
        vmin=None, vmax=None, altitude_range=(-99.9, 399.9), colorbar=False, return_image=False, annotate=True):

        if vmin is None:
            vmin = self.vmin

        if vmax is None:
            vmax = self.vmax

        if ax is None:
            ax = plt.gca()

        plt.sca(ax)
        plt.cla()
        freq_extent = (self.extent[0], self.extent[1],
            altitude_range[1], altitude_range[0])

        i = self.ionogram_list[0]
        inx = 1.0E6* (i.frequencies.shape[0] * f) / (i.frequencies[-1] - i.frequencies[0])
        img = self.tser_arr_all[:,int(inx),:]

        new_altitudes = np.arange(altitude_range[0], altitude_range[1], 14.)
        new_img = np.zeros((new_altitudes.shape[0], img.shape[1])) + np.nan

        for i in self.ionogram_list:
            e = int( round((i.time - self.extent[0]) / ais_code.ais_spacing_seconds ))

            pos = mex.iau_r_lat_lon_position(float(i.time))
            altitudes = pos[0] - ais_code.speed_of_light_kms * ais_code.ais_delays * 0.5 - mex.mars_mean_radius_km
            s = np.argsort(altitudes)
            new_img[:, e] = np.interp(new_altitudes, altitudes[s], img[s,e], left=np.nan, right=np.nan)

        plt.imshow(new_img, vmin=vmin, vmax=vmax,
            interpolation='Nearest', extent=freq_extent, origin='upper', aspect='auto')

        plt.xlim(freq_extent[0], freq_extent[1])
        plt.ylim(*altitude_range)

        ax.set_xlim(self.extent[0], self.extent[1])
        ax.xaxis.set_major_locator(celsius.SpiceetLocator())

        celsius.ylabel(r'Alt./km')
        if annotate:
            plt.annotate('f = %.1f MHz' % f, (0.02, 0.9),
                    xycoords='axes fraction', color='cyan', verticalalignment='top', fontsize='small')

        if colorbar:
            old_ax = plt.gca()
            plt.colorbar(cax = celsius.make_colorbar_cax(), ticks=self.cbar_ticks).set_label(r"$Log_{10} V^2 m^{-2} Hz^{-1}$")
            plt.sca(old_ax)

        if return_image:
            return new_img, freq_extent, new_altitudes

 def plot_lon(self, ax=None, label=True, fmt='k-', **kwargs):
     if ax is None:
         ax = plt.gca()
     plt.sca(ax)
     self.make_axis_circular(ax)
     self.generate_position()
     v = celsius.deg_unwrap(self.iau_pos[2])
     for i in [-1,0,1]:
         plt.plot(self.t, v + i * 360, fmt, **kwargs)
     if label:
         celsius.ylabel(r'$\varphi$')
     plt.ylim(0., 360.)