Python matplotlib.pyplot.fill() Examples

def make_topography_overlay_4_blockplot(self, cell_number, direction):
        p1, p2 = self.calculate_p1p2(direction, cell_number)
        resx = self.model._grid.topography.resolution[0]
        resy = self.model._grid.topography.resolution[1]
        print('p1', p1, 'p2', p2)
        x, y, z = self._slice_topo_4_sections(p1, p2, resx, resy)
        if direction == 'x':
            a = np.vstack((y, z)).T
            ext = self.model._grid.regular_grid.extent[[2, 3]]
        elif direction == 'y':
            a = np.vstack((x, z)).T
            ext = self.model._grid.regular_grid.extent[[0, 1]]
        a = np.append(a,
                      ([ext[1], a[:, 1][-1]],
                       [ext[1], self.model._grid.regular_grid.extent[5]],
                       [ext[0], self.model._grid.regular_grid.extent[5]],
                       [ext[0], a[:, 1][0]]))
        line = a.reshape(-1, 2)
        plt.fill(line[:, 0], line[:, 1], color='k') 

def plot_polygon(poly, symbol='k-', **kwargs):
    """Plots a polygon using the given symbol."""
    for i in range(poly.GetGeometryCount()):
        subgeom = poly.GetGeometryRef(i)
        x, y = zip(*subgeom.GetPoints())
        plt.plot(x, y, symbol, **kwargs)

# Use this function to fill polygons (shown shortly after
# this listing in the text). Uncomment this one and comment
# out the one above.
# def plot_polygon(poly, symbol='w', **kwargs):
#     """Plots a polygon using the given symbol."""
#     for i in range(poly.GetGeometryCount()):
#         x, y = zip(*poly.GetGeometryRef(i).GetPoints())
#         plt.fill(x, y, symbol, **kwargs)


# This function is new. 

def AreaPlot(self, Lon, Lat, Set=['y',1,'k',1]):

    x, y = self._map(Lon, Lat)

    if Set[0]:
      self._zo += 1
      plt.fill(x, y, color = Set[0],
                     alpha = Set[1],
                     zorder = self._zo)
    if Set[2]:
      self._zo += 1
      plt.plot(x, y, Set[2],
                     linewidth = Set[3],
                     zorder = self._zo)

  #--------------------------------------------------------------------------------------- 

def create_figure():
    plt.figure()
    x = np.linspace(0, 1, 15)

    # line plot
    plt.plot(x, x ** 2, "b-")

    # marker
    plt.plot(x, 1 - x**2, "g>")

    # filled paths and patterns
    plt.fill_between([0., .4], [.4, 0.], hatch='//', facecolor="lightgray",
                     edgecolor="red")
    plt.fill([3, 3, .8, .8, 3], [2, -2, -2, 0, 2], "b")

    # text and typesetting
    plt.plot([0.9], [0.5], "ro", markersize=3)
    plt.text(0.9, 0.5, u'unicode (ü, °, µ) and math ($\\mu_i = x_i^2$)',
             ha='right', fontsize=20)
    plt.ylabel('sans-serif, blue, $\\frac{\\sqrt{x}}{y^2}$..',
               family='sans-serif', color='blue')

    plt.xlim(0, 1)
    plt.ylim(0, 1)


# test compiling a figure to pdf with xelatex 

def verbosePlot(object, isPoints = False, isPaths = False, isPolygons = False):
    import numpy as np
    import matplotlib.pyplot as plt
    for child in object:
        data = np.array(child)
        if isPolygons:
            plt.fill(data.T[0], data.T[1], facecolor='grey', alpha=0.3, linestyle='--', linewidth=1)
        elif isPaths:
            plt.plot(data.T[0], data.T[1], linestyle='-', linewidth=3)
        elif isPoints:
            plt.plot(data.T[0], data.T[1], linestyle='', marker='x', markersize=10, mew=3) 

def build_mask(self):
        naxis1, naxis2 = self.wcs.pixel_shape
        edges_x = [0]*naxis2 + [naxis1-1]*naxis2 + list(range(naxis1)) * 2
        edges_y = list(range(naxis2)) * 2 + [0]*naxis1 + [naxis2-1]*naxis1

        polygon = list(zip(edges_x, edges_y))
        img = Image.new("L", (naxis1, naxis2), 0)
        ImageDraw.Draw(img).polygon(polygon, outline=1, fill=1)
        mask = np.array(img)

        self.mask = mask

#
# Utility functions used in generating or supporting the grid definitions
# 

def plot(self, output=None, color='b'):
        fig = plt.figure(figsize=(8, 6))
        ax = fig.add_subplot(111, projection="mollweide")
        ax.fill(self.corners[:, 0], self.corners[:, 1],
                facecolor='green',edgecolor='forestgreen', alpha=0.25)
        if output:
            fig.write(output)
        return ax 

def build_mask(self):
        naxis1, naxis2 = self.wcs.pixel_shape
        edges_x = [0]*naxis2 + [naxis1-1]*naxis2 + list(range(naxis1)) * 2
        edges_y = list(range(naxis2)) * 2 + [0]*naxis1 + [naxis2-1]*naxis1

        polygon = list(zip(edges_x, edges_y))
        img = Image.new("L", (naxis1, naxis2), 0)
        ImageDraw.Draw(img).polygon(polygon, outline=1, fill=1)
        mask = np.array(img)

        self.mask = mask 

def plot(self, projection='aitoff'):
        if not self.projection_cells:
            print("Please run `get_projection_cells()' first...")
            return
        plt.figure()
        plt.subplot(111, projection=projection)
        plt.grid(True)
        for pc in self.projection_cells:
            plt.fill(pc.footprint[:,0], pc.footprint[:,1],
                     facecolor='green', edgecolor='forestgreen',
                     alpha=0.25)
            plt.text(pc.footprint[0,0], pc.footprint[0,1], "{}".format(pc.cell_id),
                     horizontalalignment='right', verticalalignment='bottom') 

def plot_section_scalarfield(self, section_name, sn, levels=50, show_faults=True, show_topo=True, lithback=True):
        if self.model.solutions.sections is None:
            raise AttributeError('no sections for plotting defined')
        if self.model._grid.topography is None:
            show_topo = False
        shapes = self.model._grid.sections.resolution
        fig = plt.figure(figsize=(16, 10))
        axes = fig.add_subplot(1, 1, 1)
        j = np.where(self.model._grid.sections.names == section_name)[0][0]
        l0, l1 = self.model._grid.sections.get_section_args(section_name)
        if show_faults:
            self.extract_section_fault_lines(section_name, zorder=9)

        if show_topo:
            xy = self.make_topography_overlay_4_sections(j)
            axes.fill(xy[:, 0], xy[:, 1], 'k', zorder=10)

        axes.contour(self.model.solutions.sections[1][sn][l0:l1].reshape(shapes[j][0], shapes[j][1]).T,
                     # origin='lower',
                     levels=levels, cmap='autumn', extent=[0, self.model._grid.sections.dist[j],
                                                           self.model._grid.regular_grid.extent[4],
                                                           self.model._grid.regular_grid.extent[5]], zorder=8)
        axes.set_aspect('equal')
        if lithback:
            axes.imshow(self.model.solutions.sections[0][0][l0:l1].reshape(shapes[j][0], shapes[j][1]).T,
                        origin='lower',
                        cmap=self._cmap, norm=self._norm, extent=[0, self.model._grid.sections.dist[j],
                                                                  self.model._grid.regular_grid.extent[4],
                                                                  self.model._grid.regular_grid.extent[5]])

        labels, axname = self._make_section_xylabels(section_name, len(axes.get_xticklabels()))
        pos_list = np.linspace(0, self.model._grid.sections.dist[j], len(labels))
        axes.xaxis.set_major_locator(FixedLocator(nbins=len(labels), locs=pos_list))
        axes.xaxis.set_major_formatter(FixedFormatter((labels)))
        axes.set(title=self.model._grid.sections.names[j], xlabel=axname, ylabel='Z') 

def plot_section_by_name(self, section_name, show_data=True, show_faults=True, show_topo=True,
                             show_all_data=False, contourplot=True, radius='default', **kwargs):

        if self.model.solutions.sections is None:
            raise AttributeError('no sections for plotting defined')
        if section_name not in self.model._grid.sections.names:
            raise AttributeError(f'Section "{section_name}" is not defined. '
                                 f'Available sections for plotting: {self.model._grid.sections.names}')

        j = np.where(self.model._grid.sections.names == section_name)[0][0]
        l0, l1 = self.model._grid.sections.get_section_args(section_name)
        shape = self.model._grid.sections.resolution[j]

        image = self.model.solutions.sections[0][0][l0:l1].reshape(shape[0], shape[1]).T
        extent = [0, self.model._grid.sections.dist[j][0],
                  self.model._grid.regular_grid.extent[4], self.model._grid.regular_grid.extent[5]]

        if show_data:
            self.plot_section_data(section_name=section_name, show_all_data=show_all_data, radius=radius)

        axes = plt.gca()
        axes.imshow(image, origin='lower', zorder=-100,
                    cmap=self._cmap, norm=self._norm, extent=extent)
        if show_faults and not contourplot:
            self.extract_section_lines(section_name, axes, faults_only=True)
        else:
            self.extract_section_lines(section_name, axes, faults_only=False)
        if show_topo:
            if self.model._grid.topography is not None:
                alpha = kwargs.get('alpha', 1)
                xy = self.make_topography_overlay_4_sections(j)
                axes.fill(xy[:, 0], xy[:, 1], 'k', zorder=10, alpha=alpha)

        labels, axname = self._make_section_xylabels(section_name, len(axes.get_xticklabels()) - 1)
        pos_list = np.linspace(0, self.model._grid.sections.dist[j], len(labels))
        axes.xaxis.set_major_locator(FixedLocator(nbins=len(labels), locs=pos_list))
        axes.xaxis.set_major_formatter(FixedFormatter((labels)))
        axes.set(title=self.model._grid.sections.names[j], xlabel=axname, ylabel='Z') 

def _blob(x, y, w, w_max, area, cmap=None):
    """
    Draws a square-shaped blob with the given area (< 1) at the given coordinates.
    """
    hs = np.sqrt(area) / 2
    xcorners = np.array([x - hs, x + hs, x + hs, x - hs])
    ycorners = np.array([y - hs, y - hs, y + hs, y + hs])

    plt.fill(xcorners, ycorners, color=cmap)  # cmap(int((w + w_max) * 256 / (2 * w_max))))


# Modified from QuTip (see https://bit.ly/2LrbayH ) which in turn modified the code from the
# SciPy Cookbook. 

def create_figure():
    plt.figure()
    x = np.linspace(0, 1, 15)

    # line plot
    plt.plot(x, x ** 2, "b-")

    # marker
    plt.plot(x, 1 - x**2, "g>")

    # filled paths and patterns
    plt.fill_between([0., .4], [.4, 0.], hatch='//', facecolor="lightgray",
                     edgecolor="red")
    plt.fill([3, 3, .8, .8, 3], [2, -2, -2, 0, 2], "b")

    # text and typesetting
    plt.plot([0.9], [0.5], "ro", markersize=3)
    plt.text(0.9, 0.5, 'unicode (ü, °, µ) and math ($\\mu_i = x_i^2$)',
             ha='right', fontsize=20)
    plt.ylabel('sans-serif, blue, $\\frac{\\sqrt{x}}{y^2}$..',
               family='sans-serif', color='blue')

    plt.xlim(0, 1)
    plt.ylim(0, 1)


# test compiling a figure to pdf with xelatex 

def PlotCompTable(CompTable):

  for CT in CompTable:

    X = [CT[2], CT[3], CT[3], CT[2], CT[2]]
    Y = [CT[0], CT[0], CT[0]+CT[1], CT[0]+CT[1], CT[0]]

    plt.plot(X, Y, 'r--', linewidth=2)
    plt.fill(X, Y, color='y',alpha=0.1)

#----------------------------------------------------------------------------------------- 

def plot_optimizer(opt, lower, upper):
    import matplotlib.pyplot as plt
    plt.set_cmap("viridis")

    if not opt.models:
        print('Can not plot opt, since models do not exist yet.')
        return
    model = opt.models[-1]
    x = np.linspace(lower, upper).reshape(-1, 1)
    x_model = opt.space.transform(x)

    # Plot Model(x) + contours
    y_pred, sigma = model.predict(x_model, return_std=True)
    plt.plot(x, -y_pred, "g--", label=r"$\mu(x)$")
    plt.fill(np.concatenate([x, x[::-1]]),
             np.concatenate([-y_pred - 1.9600 * sigma,
                             (-y_pred + 1.9600 * sigma)[::-1]]),
             alpha=.2, fc="g", ec="None")

    # Plot sampled points
    plt.plot(opt.Xi, -np.array(opt.yi),
             "r.", markersize=8, label="Observations")

    # Adjust plot layout
    plt.grid()
    plt.legend(loc='best')
    plt.show() 

def negative_contour_area(mpl_obj):
    """
    Returns a array of contour levels and
    corresponding cumulative area of contours
    specifically used for calculating negative contour's area when the contours are depth of lake
    # Refer: Nikolai Shokhirev http://www.numericalexpert.com/blog/area_calculation/

    :param mpl_obj: Matplotlib contour object
    :return: [(level1, area1), (level1, area1+area2)]
    """
    n_c = len(mpl_obj.collections)  # n_c = no of contours
    print 'No. of contours = {0}'.format(n_c)
    # area = 0.0000
    cont_area_array = []
    for contour in range(n_c):
        n_p = len(mpl_obj.collections[contour].get_paths())
        zc = mpl_obj.levels[contour]
        print zc
        print n_p
        area = 0.000
        for path in range(n_p):
            p = mpl_obj.collections[contour].get_paths()[path]
            v = p.vertices
            l = len(v)
            s = 0.0000
            # plt.figure()
            # plt.fill(v[:, 0], v[:, 1], facecolor='b')
            # plt.grid()
            # plt.show()
            for i in range(l):
                j = (i + 1) % l
                s += (v[j, 0] - v[i, 0]) * (v[j, 1] + v[i, 1])
            poly_area = abs(0.5 * s)
            area += poly_area
        cont_area_array.append((zc, area))
    return cont_area_array 

def negative_contour_area(mpl_obj):
    """
    Returns a array of contour levels and
    corresponding cumulative area of contours
    specifically used for calculating negative contour's area when the contours are depth of lake
    # Refer: Nikolai Shokhirev http://www.numericalexpert.com/blog/area_calculation/

    :param mpl_obj: Matplotlib contour object
    :return: [(level1, area1), (level1, area1+area2)]
    """
    n_c = len(mpl_obj.collections)  # n_c = no of contours
    print 'No. of contours = {0}'.format(n_c)
    # area = 0.0000
    cont_area_array = []
    for contour in range(n_c):
        n_p = len(mpl_obj.collections[contour].get_paths())
        zc = mpl_obj.levels[contour]
        print zc
        print n_p
        area = 0.000
        for path in range(n_p):
            p = mpl_obj.collections[contour].get_paths()[path]
            v = p.vertices
            l = len(v)
            s = 0.0000
            # plt.figure()
            # plt.fill(v[:, 0], v[:, 1], facecolor='b')
            # plt.grid()
            # plt.show()
            for i in range(l):
                j = (i + 1) % l
                s += (v[j, 0] - v[i, 0]) * (v[j, 1] + v[i, 1])
            poly_area = abs(0.5 * s)
            area += poly_area
        cont_area_array.append((zc, area))
    return cont_area_array 

def _plot_monochromatic(observer, xy_to_2d, fill_horseshoe=True):
    # draw outline of monochromatic spectra
    lmbda = 1.0e-9 * numpy.arange(380, 701)
    values = []
    # TODO vectorize (see <https://github.com/numpy/numpy/issues/10439>)
    for k, _ in enumerate(lmbda):
        data = numpy.zeros(len(lmbda))
        data[k] = 1.0
        values.append(_xyy_from_xyz100(spectrum_to_xyz100((lmbda, data), observer))[:2])
    values = numpy.array(values)

    # Add the values between the first and the last point of the horseshoe
    t = numpy.linspace(0.0, 1.0, 101)
    connect = xy_to_2d(numpy.outer(values[0], t) + numpy.outer(values[-1], 1 - t))
    values = xy_to_2d(values.T).T
    full = numpy.concatenate([values, connect.T])

    # fill horseshoe area
    if fill_horseshoe:
        plt.fill(*full.T, color=[0.8, 0.8, 0.8], zorder=0)
    # plot horseshoe outline
    plt.plot(
        values[:, 0],
        values[:, 1],
        "-k",
        # label="monochromatic light"
    )
    # plot dotted connector
    plt.plot(connect[0], connect[1], ":k")
    return 

def plot_visible_slice(
        self, level, outline_prec=1.0e-2, plot_srgb_gamut=True, fill_color="0.8"
    ):
        # first plot the monochromatic outline
        mono_xy, conn_xy = get_mono_outline_xy(
            observer=cie_1931_2(), max_stepsize=outline_prec
        )

        mono_vals = numpy.array([self._bisect(xy, self.k0, level) for xy in mono_xy])
        conn_vals = numpy.array([self._bisect(xy, self.k0, level) for xy in conn_xy])

        k1, k2 = [k for k in [0, 1, 2] if k != self.k0]
        plt.plot(mono_vals[:, k1], mono_vals[:, k2], "-", color="k")
        plt.plot(conn_vals[:, k1], conn_vals[:, k2], ":", color="k")
        #
        if fill_color is not None:
            xyz = numpy.vstack([mono_vals, conn_vals[1:]])
            plt.fill(xyz[:, k1], xyz[:, k2], facecolor=fill_color, zorder=0)

        if plot_srgb_gamut:
            self._plot_srgb_gamut(self.k0, level)

        plt.axis("equal")
        plt.xlabel(self.labels[k1])
        plt.ylabel(self.labels[k2])
        plt.title(f"{self.labels[self.k0]} = {level}") 

def _plot_sequences(self, seq_ind, seqs, ro_inds, x_unit, bounds):
        """
        The actual plotting of the sequences happens here.
        
        Args:
            seqs:            list of sequences to be plotted (i.e. one list of sequence for each channel)
            ro_inds:         indices of the readout
            x_unit:          x_unit for the time axis
            bounds:          boundaries of the plot (xmin, xmax, ymin, ymax)
            show_quadrature: set to "I" or "Q" if you want to display either quadrature instead of the amplitude
        """
        if not qkit.module_available("matplotlib"):
            raise ImportError("matplotlib not found.")

        fig = plt.figure(figsize=(18,6))
        xmin, xmax, ymin, ymax = bounds
        samplerate = self._sample.clock
        # plot sequence
        for i, chan in enumerate(self.channels[1:]):
            if len(ro_inds[i]) > seq_ind:
                chan._plot_sequence(seqs[i][seq_ind], ro_inds[i][seq_ind], x_unit, bounds, col = self._chancols[i + 1], plot_readout = False)
        
            # plot readout
        plt.fill([0, 0, - self._sample.readout_tone_length / self._x_unit[x_unit], - self._sample.readout_tone_length / self._x_unit[x_unit]], 
                [0, ymax, ymax, 0], color = "C7", alpha = 0.3)
        # add label for readout
        plt.text(-0.5*self._sample.readout_tone_length / self._x_unit[x_unit], ymax/2.,
                "readout", horizontalalignment = "center", verticalalignment = "center", rotation = 90, size = 14)
        
        # adjust plot limits
        plt.xlim(xmin + 0.005 * abs(xmax - xmin), xmax - 0.006 * abs(xmax - xmin))
        plt.ylim(ymin, ymax + 0.025 * (ymax - ymin))
        plt.xlabel("time " + x_unit)
        return 

def cl_dist_map(x,y,z,xmin,xmax,ymin,ymax,dx):
    """function for centerline rasterization and distance map calculation (does not return zmap)
    used for cutoffs only 
    inputs:
    x,y,z - coordinates of centerline
    xmin, xmax, ymin, ymax - x and y coordinates that define the area of interest
    dx - gridcell size (m)
    returns:
    cl_dist - distance map (distance from centerline)
    x_pix, y_pix, - x and y pixel coordinates of the centerline
    """
    y = y[(x>xmin) & (x<xmax)]
    z = z[(x>xmin) & (x<xmax)]
    x = x[(x>xmin) & (x<xmax)]    
    xdist = xmax - xmin
    ydist = ymax - ymin
    iwidth = int((xmax-xmin)/dx)
    iheight = int((ymax-ymin)/dx)
    xratio = iwidth/xdist
    # create list with pixel coordinates:
    pixels = []
    for i in range(0,len(x)):
        px = int(iwidth - (xmax - x[i]) * xratio)
        py = int(iheight - (ymax - y[i]) * xratio)
        pixels.append((px,py))
    # create image and numpy array:
    img = Image.new("RGB", (iwidth, iheight), "white")
    draw = ImageDraw.Draw(img)
    draw.line(pixels, fill="rgb(0, 0, 0)") # draw centerline as black line
    pix = np.array(img)
    cl = pix[:,:,0]
    cl[cl==255] = 1 # set background to 1 (centerline is 0)
    # calculate Euclidean distance map:
    cl_dist, inds = ndimage.distance_transform_edt(cl, return_indices=True)
    y_pix,x_pix = np.where(cl==0)
    return cl_dist, x_pix, y_pix 

def _plot_sequence(self, seq, ro_ind, x_unit, bounds, col = "C0", plot_readout = True):
        """
        The actual plotting of the sequences happens here.
        
        Input:
            seq:    sequence (as array) to be plotted
            ro_ind: index of the readout in seq
            x_unit: x_unit for the time axis
            bounds: boundaries for the plot (xmin, xmax, ymin, ymax)
        """
        if not qkit.module_available("matplotlib"):
            raise ImportError("matplotlib not found.")
        
        if plot_readout:
            fig = plt.figure(figsize = (18, 6))
        xmin, xmax, ymin, ymax = bounds
        samplerate = self._sample.clock
        time = -(np.arange(0, len(seq) + 1, 1) - ro_ind) / (samplerate * self._x_unit[x_unit])
        # make sure last point of the waveform goes to zero
        seq = np.append(seq, 0)

        # plot sequence
        plt.plot(time, seq, col)
        plt.fill(time, seq, color = col, alpha = 0.2)
        # plot readout
        if plot_readout:
            plt.fill([0, 0, - self._sample.readout_tone_length / self._x_unit[x_unit], - self._sample.readout_tone_length / self._x_unit[x_unit]], 
                    [0, ymax, ymax, 0], color = "C7", alpha = 0.3)        
            # add label for readout
            plt.text(-0.5*self._sample.readout_tone_length / self._x_unit[x_unit], ymax/2.,
                    "readout", horizontalalignment = "center", verticalalignment = "center", rotation = 90, size = 14)
        # adjust bounds
        plt.xlim(xmin + 0.005 * abs(xmax - xmin), xmax - 0.006 * abs(xmax - xmin))
        plt.ylim(ymin, ymax + 0.025 * (ymax - ymin))
        plt.xlabel("time " + x_unit)
        return 

def display(x, color, list_save=None) :
    kde  = KernelDensity(kernel='gaussian', bandwidth= .005 ).fit(x.data.cpu().numpy())
    dens = np.exp( kde.score_samples(t_plot) )
    dens[0] = 0 ; dens[-1] = 0;
    plt.fill(t_plot, dens, color=color)
    if list_save is not None :
        list_save.append(dens.ravel()) # We'll save a csv at the end 

def create_figure():
    plt.figure()
    x = np.linspace(0, 1, 15)

    # line plot
    plt.plot(x, x ** 2, "b-")

    # marker
    plt.plot(x, 1 - x**2, "g>")

    # filled paths and patterns
    plt.fill_between([0., .4], [.4, 0.], hatch='//', facecolor="lightgray",
                     edgecolor="red")
    plt.fill([3, 3, .8, .8, 3], [2, -2, -2, 0, 2], "b")

    # text and typesetting
    plt.plot([0.9], [0.5], "ro", markersize=3)
    plt.text(0.9, 0.5, 'unicode (ü, °, µ) and math ($\\mu_i = x_i^2$)',
             ha='right', fontsize=20)
    plt.ylabel('sans-serif, blue, $\\frac{\\sqrt{x}}{y^2}$..',
               family='sans-serif', color='blue')

    plt.xlim(0, 1)
    plt.ylim(0, 1)


# test compiling a figure to pdf with xelatex 

def poly_area(xy):
    """
    Calculates polygon area
    x = xy[:,0], y[xy[:,1]
    :param xy:
    :return:
    """
    l = len(xy)
    s = 0.0
    for i in range(l):
        j = (i+1) % l
        s += (xy[j, 0] - xy[i, 0]) * (xy[j,1]+ xy[i,1])
    return -0.5*s

# # zero contour has two paths 0, 1
# p_0_0 = CS.collections[0].get_paths()[0]    # CS.collections[index of contour].get_paths()[index of path]
# p_0_1 = CS.collections[0].get_paths()[1]
# v_0_0 = p_0_0.vertices
# v_0_1 = p_0_1.vertices
# area_0_0 = abs(poly_area(v_0_0))
# area_0_1 = abs(poly_area(v_0_1))
# area_0 = area_0_0 + area_0_1
# z_0 = CS.levels[0]
# # print z_0, area_0
# plt.fill(v_0_0[:,0], v_0_0[:,1], facecolor='g')
# plt.show()
# # 0.4 contour has three paths 0,1,2
# print len(CS.collections[21].get_paths()) 

def plot_all_sections(self, show_data=False, section_names=None, show_topo=True,
                          figsize=(12, 12)):
        if self.model.solutions.sections is None:
            raise AttributeError('no sections for plotting defined')
        if self.model._grid.topography is None:
            show_topo = False
        if section_names is not None:
            if isinstance(section_names, list):
                section_names = np.array(section_names)
        else:
            section_names = self.model._grid.sections.names

        shapes = self.model._grid.sections.resolution
        fig, axes = plt.subplots(nrows=len(section_names), ncols=1, figsize=figsize)
        for i, section in enumerate(section_names):
            j = np.where(self.model._grid.sections.names == section)[0][0]
            l0, l1 = self.model._grid.sections.get_section_args(section)

            self.extract_section_lines(section, axes[i], faults_only=False)

            if show_topo:
                xy = self.make_topography_overlay_4_sections(j)
                axes[i].fill(xy[:, 0], xy[:, 1], 'k', zorder=10)

            # if show_data:
            #    section = str(section)
            #    print(section)
            #    self.plot_section_data(section_name=section)

            axes[i].imshow(self.model.solutions.sections[0][0][l0:l1].reshape(shapes[j][0], shapes[j][1]).T,
                           origin='lower', zorder=-100,
                           cmap=self._cmap, norm=self._norm, extent=[0, self.model._grid.sections.dist[j],
                                                                     self.model._grid.regular_grid.extent[4],
                                                                     self.model._grid.regular_grid.extent[5]])

            labels, axname = self._make_section_xylabels(section, len(axes[i].get_xticklabels()) - 1)
            pos_list = np.linspace(0, self.model._grid.sections.dist[j], len(labels))
            axes[i].xaxis.set_major_locator(FixedLocator(nbins=len(labels), locs=pos_list))
            axes[i].xaxis.set_major_formatter(FixedFormatter((labels)))
            axes[i].set(title=self.model._grid.sections.names[j], xlabel=axname, ylabel='Z')

        fig.tight_layout() 

def EStarRStarResiduals(Sc=0.71):

    """
    MDH
    
    """

    # setup the figure
    Fig = CreateFigure(FigSizeFormat="EPSL",AspectRatio=4.)
    Ax = plt.subplot(111)
    
    # Calculate analytical relationship
    EStar = np.logspace(-1,3,1000)
    RStar = CalculateRStar(EStar)
    
    plt.plot([-6,41],[0,0],'k--',lw=0.5,zorder=1)
    plt.fill([-6,41,41,-6],[0,0,0.5,0.5],color=[1.0,0.8,0.8],zorder=0)
    plt.fill([-6,41,41,-6],[0,0,-0.5,-0.5],color=[0.8,0.9,1.0],zorder=0)
    plt.text(-5,0.35,"Growing")
    plt.text(-5,-0.45,"Decaying")
    
    #loop through the basins
    for Basin in range(0,NoBasins):
    #for Basin in range(0,1):

        # Get the hillslope data for the basin        
        Data = CalculateEStarRStar(Basin)
        
        # Calculate log transformed orthogonal residuals
        Residuals, Xortho, Yortho = OrthogonalResiduals(np.log10(EStar), np.log10(RStar), np.log10(Data.EStar.as_matrix()), np.log10(Data.RStar.as_matrix()))

        # plot violin of residuals
        violin_parts = plt.violinplot(-Residuals, [Basin,],showmeans=False, showmedians=False, showextrema=False)
        
        # set the colour of the distribution
        for fc in violin_parts['bodies']:
            fc.set_facecolor([0.2,0.2,0.2])
        
        # plot the median
        Median = -np.median(Residuals)
        plt.plot([Basin-0.2,Basin+0.2],[Median,Median],'k-',lw=1)
            
        
        
    plt.xlabel("Basin Number")
    plt.ylabel("Orthogonal Residuals\n($log_{10}\:E* R*$)")
    plt.xlim(-6,41)
    plt.ylim(-0.5,0.5)
    plt.xticks(np.arange(0,NoBasins,2))
    Ax.tick_params(axis='x',which='minor',bottom='on')
    
    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.savefig(PlotDirectory+FilenamePrefix + "_ESRSOrthoResidualsViolin.png", dpi=300) 

def polygon(polygon, style='-k', linewidth=1, fill=None, alpha=1., label=None,
            xy2ne=False, linealpha=1.):
    """
    Plot a polygon.

    Parameters:

    * polygon : :class:`geoist.inversion.geometry.Polygon`
        The polygon
    * style : str
        Color and line style string (as in matplotlib.pyplot.plot)
    * linewidth : float
        Line width
    * fill : str
        A color string used to fill the polygon. If None, the polygon is not
        filled
    * alpha : float
        Transparency of the fill (1 >= alpha >= 0). 0 is transparent and 1 is
        opaque
    * linealpha : float
        Transparency of the line (1 >= alpha >= 0). 0 is transparent and 1 is
        opaque
    * label : str
        String with the label identifying the polygon in the legend
    * xy2ne : True or False
        If True, will exchange the x and y axis so that the x coordinates of
        the polygon are north. Use this when drawing on a map viewed from
        above. If the y-axis of the plot is supposed to be z (depth), then use
        ``xy2ne=False``.

    Returns:

    * lines : matplotlib Line object
        Line corresponding to the polygon plotted

    """
    if xy2ne:
        tmpx = [y for y in polygon.y]
        tmpx.append(polygon.y[0])
        tmpy = [x for x in polygon.x]
        tmpy.append(polygon.x[0])
    else:
        tmpx = [x for x in polygon.x]
        tmpx.append(polygon.x[0])
        tmpy = [y for y in polygon.y]
        tmpy.append(polygon.y[0])
    kwargs = {'linewidth': linewidth, 'alpha': linealpha}
    if label is not None:
        kwargs['label'] = label
    line, = pyplot.plot(tmpx, tmpy, style, **kwargs)
    if fill is not None:
        pyplot.fill(tmpx, tmpy, color=fill, alpha=alpha)
    return line 

def square(area, style='-k', linewidth=1, fill=None, alpha=1., label=None,
           xy2ne=False):
    """
    Plot a square.

    Parameters:

    * area : list = [x1, x2, y1, y2]
        Borders of the square
    * style : str
        String with the color and line style (as in matplotlib.pyplot.plot)
    * linewidth : float
        Line width
    * fill : str
        A color string used to fill the square. If None, the square is not
        filled
    * alpha : float
        Transparency of the fill (1 >= alpha >= 0). 0 is transparent and 1 is
        opaque
    * label : str
        label associated with the square.
    * xy2ne : True or False
        If True, will exchange the x and y axis so that the x coordinates of
        the polygon are north. Use this when drawing on a map viewed from
        above. If the y-axis of the plot is supposed to be z (depth), then use
        ``xy2ne=False``.

    Returns:

    * axes : ``matplitlib.axes``
        The axes element of the plot

    """
    x1, x2, y1, y2 = area
    if xy2ne:
        x1, x2, y1, y2 = y1, y2, x1, x2
    xs = [x1, x1, x2, x2, x1]
    ys = [y1, y2, y2, y1, y1]
    kwargs = {'linewidth': linewidth}
    if label is not None:
        kwargs['label'] = label
    plot, = pyplot.plot(xs, ys, style, **kwargs)
    if fill is not None:
        pyplot.fill(xs, ys, color=fill, alpha=alpha)
    return plot 

def fill_between(xdata, lb, ub, xlabel='', ylabel='', alpha=0.3, color=None):
    """Fill between lines.

    Parameters
    ----------
    xdata : |ndarray|
        X values of points to plot.  Should be vector of length ``Nsamples``.
    lb : |ndarray|
        Lower bound of fill.  Should be of size ``(Nsamples,...)``.
    ub : |ndarray|
        Upper bound of fill.  Should be same size as lb.
    xlabel : str
        Label for the x axis. Default is no x axis label.
    ylabel : str
        Label for the y axis.  Default is no y axis label.
    fmt : str or matplotlib linespec
        Line marker to use.  Default = ``'-'`` (a normal line).
    color : matplotlib color code or list of them
        Color(s) to use to plot the distribution.
        See https://matplotlib.org/tutorials/colors/colors.html
        Default = use the default matplotlib color cycle
    """

    # Check shapes
    if not np.all(lb.shape == ub.shape):
        raise ValueError('lb and ub must have same shape')
    if len(xdata) != lb.shape[0]:
        raise ValueError('xdata does not match shape of lb and ub')

    # If 1d make 2d
    if lb.ndim == 1:
        lb = np.expand_dims(lb, 1)
        ub = np.expand_dims(ub, 1)

    # Number of fills and datasets
    dims = lb.shape[1:]
    Nd = int(np.prod(dims))
    Np = lb.shape[0]

    # Flatten if >1D
    lb = np.reshape(lb, (lb.shape[0], Nd), order='F')
    ub = np.reshape(ub, (ub.shape[0], Nd), order='F')

    # Plot the data
    for iD in range(Nd): #for each dataset,
        next_color = get_next_color(color, iD)
        lab = get_ix_label(iD, dims)
        plt.fill_between(xdata, lb[:,iD], ub[:,iD],
                         alpha=alpha, facecolor=next_color,
                         label=lab)

    # Only show the legend if there are >1 datasets
    if Nd > 1:
        plt.legend()

    # Set x axis label, and no y axis or bounding box needed
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)