ordplot

一个gps小工具包

    xMajorFmt=pylab.DateFormatter("%H:%M\n%j")
    xMinorFmt=pylab.NullFormatter()
    xMajorLoc=matplotlib.dates.DayLocator()
    xMinorLoc=matplotlib.dates.HourLocator()

    rExtent=0.89
    if options.legend:
        rExtent=0.82

    if axesCount == 2:
        ax1 = fig.add_axes([0.08, 0.52, rExtent, 0.42])
    elif axesCount == 1:
        ax1 = fig.add_axes([0.08, 0.10, rExtent, 0.82])
    ax1.set_ylim((0,0))

    if plotOrds:
        label="All prns"
        for prn in prns:
            onePrn = ords[:,ords[1]==prn]
            oneWart = warts[:,warts[1]==prn]
            if prn in highlightPrns:
                smoothPlot(ax1, onePrn[0], onePrn[2], ',', "prn %2d"%prn,
                           options.smoothing)
                if options.wartiness:
                    smoothPlot(ax1, oneWart[0], oneWart[2], ',', "prn %2dW"%prn, 0)
            else:
                smoothPlot(ax1, ords[0], ords[2], 'g,', label, 0)
                if options.wartiness:
                    smoothPlot(ax1, oneWart[0], oneWart[2], 'r,', label, 0)
                label="none"

        setLimit(ax1, ords[2], options.yRange)

        if options.clockDelta:
            # Plot the first deritive of the clock...
            ax3=fig.add_axes(ax1.get_position())
            time_diff = 86400 * (clocks[0,1:] - clocks[0,0:-1])
            clock_diff = 1e9 * 3.3356e-9 * (clocks[1,1:] - clocks[1,0:-1])
            dc = 1e3 * clock_diff/time_diff
            smoothPlot(ax1, clocks[0,:-1], dc, 'b-', "clk dif"%prn,
                       options.smoothing, options.yRange)
            setLimit(ax1, dc, options.yRange)
        elif len(ocds[0]):
            # otherwise plot the clock deviations
            ax3=fig.add_axes(ax1.get_position())
            smoothPlot(ax3, ocds[0], ocds[1], 'b-', "clk res"%prn,
                       options.smoothing, options.yRange)
            setLimit(ax3, ocds[1], options.yRange)
            
        # If there are rle clocks, draw a vertical line where each new model
        # starts
        if not options.clockDelta:
            for t0, t1, y0, m, d in rleClockList:
                ax1.axvline(t0, label='_nolegend_')

        if options.legend:
            ax1.legend(numpoints=2, pad=0.1, labelsep = 0, handlelen=0.005,
                       handletextsep=0.01, axespad=0.0, loc=(1,0))
            leg = pylab.gca().get_legend()
            ltext = leg.get_texts()
            llines = leg.get_lines()
            lframe = leg.get_frame()
            lframe.set_facecolor('0.4')
            pylab.setp(ltext, size=8, family="sans-serif")
            pylab.setp(llines, linewidth=2)
            leg.draw_frame(False)
        if not options.clockDelta:
            ax1.set_ylabel('ord (meters)', **yprops)
        else:
            ax1.set_ylabel('ord (meters)\nclock diff (10^-12)', **yprops)
    
        ax1.grid(True)
        pylab.figtext(rExtent+.08, 0.95, "y range %s"%options.yRange,
                      **scale_props)

    # This allows the creation of futher axes that will share the x axis
    # with the first plot.
    axprops = dict()
    axprops['sharex'] = ax1

    if axesCount == 2:
        ax2 = fig.add_axes([0.08, 0.10, rExtent, 0.38], **axprops)
    elif axesCount == 1:
        ax2 = ax1

    if plotClocks:
        if options.debugLevel:
            print "Debiasing clocks"
        clockOffset = int(numpy.average(clocks[1]))
        clocks[1] = clocks[1] - clockOffset

        pylab.figtext(rExtent+0.08, 0.49, 
                      "Clock offset removed: %d m"%clockOffset,
                      **scale_props)

        smoothPlot(ax2, clocks[0], clocks[1], 'g,', "clock offset",
                   options.smoothing, options.yRange)
        ax2.grid(True)
        ax2.set_ylabel('clock (meters)', **yprops)

        # Only plot the linear clock estimate if there is data for it...
        for t0, t1, y0, m, d in rleClockList:
            y1 = y0 + m * (t1 - t0)
            t = numpy.array([t0, t1])
            y = numpy.array([y0, y1])
            ax2.plot_date(t, y, 'b-', linewidth=1, label='_nolegend_')
            yu = y + d
            yl = y - d
            yy = pylab.concatenate( (yu, yl[::-1]) )
            tt = pylab.concatenate( (t, t[::-1]) )
            ax2.fill(tt, yy, facecolor='b', alpha=0.4, label='_nolegend_')
            
        if options.legend:
            ax2.legend(numpoints=2, pad=0.1, labelsep = 0, handlelen=0.005,
                       handletextsep=0.01, axespad=0.0, loc=(1,0))
            leg = pylab.gca().get_legend()
            pylab.setp(leg.get_texts(), size=8, family="sans-serif")
            leg.draw_frame(False)


    ax2.xaxis.set_major_formatter(xMajorFmt)
    xlabels=ax2.get_xticklabels()
    ylabels=ax2.get_yticklabels()
    pylab.setp(xlabels, fontsize=10, family='sans-serif')
    pylab.setp(ylabels, fontsize=10, family='sans-serif')

    ax2.xaxis.set_minor_formatter(xMinorFmt)
    
    # set x axis range
    if options.tStart:
        tMin = parse_time(options.tStart.split()[0:3])
    else:
        if len(ords[0]):
            tMin = min(ords[0])
        else:
            tMin = min(clocks[0])
            
    if options.tEnd:
        tMax = parse_time(options.tEnd.split()[0:3])
    else:
        if len(ords[0]):
            tMax = max(ords[0])
        else:
            tMax = max(clocks[0])
        
    ax2.set_xlim(xmin=tMin, xmax=tMax)

    # Axis labels on the upper plot would be bad since they would be
    # drawn over the upper part of the lower plot
    if axesCount > 1:
        pylab.setp(ax1.get_xticklabels(), visible=False)
    ax1.set_title(options.title)

    if (options.saveFig == None):
        pylab.show()
    else:
       pylab.savefig(options.saveFig)
# end of main


def parse_time(words):
    fsec = float(words[2][8:10])
    ydhms =  words[0]+" "+words[1]+" "+words[2][0:8]
    utime = time.strptime(ydhms, "%Y %j %H:%M:%S")
    dtime = datetime.datetime(utime[0], utime[1], utime[2],
                              utime[3], utime[4], utime[5], int(fsec*1e6))
    t0 = matplotlib.dates.date2num(dtime)
    return t0
# end of parse_time()


def setLimit(ax, y, limit=None):
    if not limit:
        lim = (min(y), max(y))
    elif limit == "ignore":
        lim = ax.get_ylim()
    elif limit[-1] == "%":
        fraction = float(limit[:-1])/100.0
        drop = int(numpy.floor(0.5 * (1-fraction) * len(y)))
        ds = numpy.sort(y)
        lim = (ds[drop], ds[-drop])
    else:
        lim = (-float(limit), float(limit))
    ax.set_ylim(lim)


def smoothPlot(ax, x, y, style="", txt="", window_len=0, limit=None):
    if x.ndim != 1 or y.ndim !=1:
        raise ValueError, "Only accepts rank 1 arrays."

    if x.size<2 or y.size<2:
        return

    n = window_len/(86400 * (x[1] - x[0]))
    if n>3 and n<y.size:
        s=numpy.r_[2*y[0]-y[n:1:-1],y,2*y[-1]-y[-1:-n:-1]]
        w=eval('numpy.bartlett(n)')
        y=numpy.convolve(w/w.sum(),s,mode='same')
        y=y[n-1:-n+1]
        y=y[n-1:-n+1]
        x=x[n-1:-n+1]

    if txt == "none":
        ax.plot_date(x, y, style)
    else:
        ax.plot_date(x, y, style, label=txt)


if __name__ == "__main__":
    main()