def write_body(self):
        try:
            from matplotlib.path import Path
        except ImportError:
            Path = None
            from matplotlib.patches import Circle, Polygon
        else:
            from matplotlib.patches import Circle, PathPatch

        indices = self.X[:, 2].argsort()
        for a in indices:
            xy = self.X[a, :2]
            if a < self.natoms:
                r = self.d[a] / 2
                if ((xy[1] + r > 0) and (xy[1] - r < self.h) and
                    (xy[0] + r > 0) and (xy[0] - r < self.w)):
                    circle = Circle(xy, r, facecolor=self.colors[a])
                    circle.draw(self.renderer)
            else:
                a -= self.natoms
                c = self.T[a]
                if c != -1:
                    hxy = self.D[c]
                    if Path is None:
                        line = Polygon((xy + hxy, xy - hxy))
                    else:
                        line = PathPatch(Path((xy + hxy, xy - hxy)))
                    line.draw(self.renderer)

 def __init__(self, xyz, **kargs):
     self._gl_3dpath = kargs.pop('gl_3dpath', None)
     self._gl_lighting = kargs.pop('gl_lighting', True)
     xy = xyz[:,0:2]
     Polygon.__init__(self, xy, **kargs)
     self.do_stencil_test = True
     self.set_3d_properties(zs = xyz[:,2])

    def __init__(self, x, y, dx, dy, width=0.001, length_includes_head=False, \
        head_width=None, head_length=None, shape='full', overhang=0, \
        head_starts_at_zero=False,**kwargs):
        """Returns a new Arrow.

        length_includes_head: True if head is counted in calculating the length.
        
        shape: ['full', 'left', 'right']
        
        overhang: distance that the arrow is swept back (0 overhang means
        triangular shape).

        head_starts_at_zero: if True, the head starts being drawn at coordinate
        0 instead of ending at coordinate 0.
        """
        if head_width is None:
            head_width = 3 * width
        if head_length is None:
            head_length = 1.5 * head_width

        distance = sqrt(dx**2 + dy**2)
        if length_includes_head:
            length=distance
        else:
            length=distance+head_length
        if not length:
            verts = [] #display nothing if empty
        else:
            #start by drawing horizontal arrow, point at (0,0)
            hw, hl, hs, lw = head_width, head_length, overhang, width
            left_half_arrow = array([
                [0.0,0.0],                  #tip
                [-hl, -hw/2.0],             #leftmost
                [-hl*(1-hs), -lw/2.0], #meets stem
                [-length, -lw/2.0],          #bottom left
                [-length, 0],
            ])
            #if we're not including the head, shift up by head length
            if not length_includes_head:
                left_half_arrow += [head_length, 0]
            #if the head starts at 0, shift up by another head length
            if head_starts_at_zero:
                left_half_arrow += [head_length/2.0, 0]
            #figure out the shape, and complete accordingly
            if shape == 'left':
                coords = left_half_arrow
            else:
                right_half_arrow = left_half_arrow*[1,-1]
                if shape == 'right':
                    coords = right_half_arrow
                elif shape == 'full':
                    coords=concatenate([left_half_arrow,right_half_arrow[::-1]])
                else:
                    raise ValueError, "Got unknown shape: %s" % shape
            cx = float(dx)/distance
            sx = float(dy)/distance
            M = array([[cx, sx],[-sx,cx]])
            verts = dot(coords, M) + (x+dx, y+dy)
        
        Polygon.__init__(self, map(tuple, verts), **kwargs)

 def __init__(self, xy, closed=True, plotview=None, **opts):
     shape = Polygon(xy, closed, **opts)
     if 'linewidth' not in opts:
         shape.set_linewidth(1.0)
     if 'facecolor' not in opts:
         shape.set_facecolor('r')
     super(NXpolygon, self).__init__(shape, resize=False, plotview=plotview)
     self.shape.set_label('Polygon')
     self.polygon = self.shape

def getGuards(polygon_no):
	positions = []
	data = triangulate(polygons[str(polygon_no)])
	vertices = data['tri']['vertices']
	colours = get_colour(data)
	min_c = get_min_colour(colours)
	stars = get_stars(data)
	star_polygons = star_poly_array(stars, vertices)
	# print(stars)
	# print(star_polygons)
	# print(vertices)
	# print(min_c)
	adj = adj_matrix(data)
	# print(adj[24])
	# print(adj[33])
	# print(adj[60])
	# pdb.set_trace()

	def compare(plt, A, B): 
	    ax1 = plt.subplot(121, aspect='equal')
	    triangle.plot.plot(ax1, **A)
	    ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
	    triangle.plot.plot(ax2, **B)
	    return (ax1, ax2)

	ax1, ax2 = compare(plt, data['original'], data['tri'])

	# Colour polygons
	patches = []
	for polygon in star_polygons:
		p = Polygon(polygon, True)
		patches.append(p)
	p = PatchCollection(patches, cmap=matplotlib.cm.jet)
	plt.gca().add_collection(p)

	# Label vertex indices
	for i, v in enumerate(vertices):
		ax1.text(v[0], v[1], str(i))

	for c, v in zip(colours, vertices):
		if c == min_c:
			positions.append(v.tolist())
		plt.annotate(c, xy=tuple(v), color='darkblue')

	# positions = starise.kernels(stars, vertices)

	# Plot kernel points
	ax1.plot(*zip(*positions), marker='^', color='g', ls='')

	# Set colours and config for plot
	colors = range(0,1000, int(1000/len(patches)))
	p.set_array(np.array(colors))	
	plt.subplots_adjust(left=0, bottom=0, right=1, top=1,
                wspace=0, hspace=0)
	plt.show()
	# pdb.set_trace()
	return positions

 def plotReachSet_norm1(self, NUM, figname):
     fig = p.figure()
     for j in range(n):
         ax = fig.add_subplot(2,2,j+1 , aspect='equal')
         ax.set_xlim(0, 4)
         ax.set_ylim(0, 1)
         ax.set_xlabel('$x_'+str(j+1)+'$')
         ax.set_ylabel('$y_'+str(j+1)+'$')
         for trace in self:
             for i in [int(floor(k*len(trace.T)/NUM)) for k in range(NUM)]:
                 verts = [(trace.x[i][j] + 1/trace.d_norm1[i][2*j], trace.y[i][j]                          ),
                          (trace.x[i][j]                            , trace.y[i][j] - 1/trace.d_norm1[i][2*j+1]),
                          (trace.x[i][j] - 1/trace.d_norm1[i][2*j], trace.y[i][j]                          ),
                          (trace.x[i][j]                              , trace.y[i][j] + 1/trace.d_norm1[i][2*j+1])]
                 # poly = Ellipse((trace.x[i][j],trace.y[i][j]), width=trace.d1[i], height=trace.d2[i], angle=trace.theta[i])
                 poly = Polygon(verts, facecolor='0.8', edgecolor='k')
                 ax.add_artist(poly)
                 poly.set_clip_box(ax.bbox)
                 poly.set_alpha(1)
                 if i==0:
                     poly.set_facecolor('r')
                 else:
                     poly.set_facecolor(p.rand(3))
         #for trace in self:
             #e = Ellipse((trace.x[0][j],trace.y[0][j]), width=trace.d1[0], height=trace.d2[0], angle=trace.theta[0])
             #ax.add_artist(e)
             #e.set_clip_box(ax.bbox)
             #e.set_alpha(1)
             #e.set_facecolor('r')
 #e.set_edgecolor('r')
     p.savefig(figname)

def get_features_from_shape(basemap, path = 'data/shape/contour_bv_MRCC/Bassins_MRCC_utm18.shp', linewidth = 2,
                        edgecolor = 'k', face_color = "none", id_list = None, zorder = 0, alpha = 1):
    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataStore = driver.Open(path, 0)
    layer = dataStore.GetLayer(0)
    latlong = osr.SpatialReference()
    latlong.ImportFromProj4("+proj=latlong")
    result = []

    id_list_lower = []
    if id_list is not None:
        id_list_lower = map(lambda x: x.lower(), id_list)

    feature = layer.GetNextFeature()
    while feature:
        geom = feature.GetGeometryRef()
        geom.TransformTo(latlong)


        #get fields of the feature
#        for i in xrange(feature.GetFieldCount()):
#            print feature.GetField(i)

#        print geom.ExportToWkt()

        polygon = loads(geom.ExportToWkt())
        boundary = polygon.exterior
        coords = np.zeros(( len(boundary.coords), 2))
        currentId = None
        for i, the_coord in enumerate(boundary.coords):


#            if feature.GetFieldAsString('abr').lower() == 'rdo':
#                print the_coord[0], the_coord[1]
            if basemap is not None:
                coords[i, 0], coords[i, 1] = basemap( the_coord[0], the_coord[1] )

            currentId = feature.GetFieldAsString("abr").lower()

        to_add = True
        if id_list is not None:
            to_add = currentId in id_list_lower

        if to_add:
            p = Polygon(coords,linewidth = linewidth, edgecolor = edgecolor,
                facecolor=face_color, zorder = zorder, alpha = alpha)
            p.basin_id = currentId
            result.append(p)
        feature = layer.GetNextFeature()


    dataStore.Destroy()
    return result

 def add_to_axes(self, ax=None, **kwargs):
     polys = [(c, self.fields[c]) for c in self.fclasses if self.fields[c]['poly']]
     if ax is None:
         fig, ax = plt.subplots(1)
     pgns = []
     for c, f in polys:
         label = f['names']
         pg = Polygon(f['poly'], closed=True, **kwargs)
         pgns.append(pg)
         x, y = pg.get_xy().T
         ax.annotate(c, xy=(np.nanmean(x), np.nanmean(y)))
         ax.add_patch(pg)

    def __init__(self, ax, poly_xy=None, max_ds=10):
        self.showverts = True
        self.max_ds = max_ds
        if poly_xy is None:
            poly_xy = default_vertices(ax)
        self.poly = Polygon(poly_xy, animated=True,
                            fc='y', ec='none', alpha=0.4)

        ax.add_patch(self.poly)
        ax.set_clip_on(False)
        ax.set_title("Click and drag a point to move it; "
                     "'i' to insert; 'd' to delete.\n"
                     "Close figure when done.")
        self.ax = ax

        x, y = zip(*self.poly.xy)
        self.line = plt.Line2D(x, y, color='none', marker='o', mfc='r',
                               alpha=0.2, animated=True)
        self._update_line()
        self.ax.add_line(self.line)

        self.poly.add_callback(self.poly_changed)
        self._ind = None # the active vert

        canvas = self.poly.figure.canvas
        canvas.mpl_connect('draw_event', self.draw_callback)
        canvas.mpl_connect('button_press_event', self.button_press_callback)
        canvas.mpl_connect('button_release_event', self.button_release_callback)
        canvas.mpl_connect('key_press_event', self.key_press_callback)
        canvas.mpl_connect('motion_notify_event', self.motion_notify_callback)
        self.canvas = canvas

    def __init__(self, ax, poly_xy=None, max_ds=10, line_width=2,
                 line_color=(0, 0, 1), face_color=(1, .5, 0)):
        self.showverts = True
        self.max_ds = max_ds
        if poly_xy is None:
            poly_xy = default_vertices(ax)
        self.poly = Polygon(poly_xy, animated=True,
                            fc=face_color, ec='none', alpha=0.4)

        ax.add_patch(self.poly)
        ax.set_clip_on(False)
        ax.set_title("Click and drag a point to move it; "
                     "'i' to insert; 'd' to delete.\n"
                     "Close figure when done.")
        self.ax = ax

        x, y = zip(*self.poly.xy)
        #line_color = 'none'
        color = np.array(line_color) * .6
        marker_face_color = line_color
        line_kwargs = {'lw': line_width, 'color': color, 'mfc': marker_face_color}
        self.line = plt.Line2D(x, y, marker='o', alpha=0.8, animated=True, **line_kwargs)
        self._update_line()
        self.ax.add_line(self.line)

        self.poly.add_callback(self.poly_changed)
        self._ind = None  # the active vert

        canvas = self.poly.figure.canvas
        canvas.mpl_connect('draw_event', self.draw_callback)
        canvas.mpl_connect('button_press_event', self.button_press_callback)
        canvas.mpl_connect('button_release_event', self.button_release_callback)
        canvas.mpl_connect('key_press_event', self.key_press_callback)
        canvas.mpl_connect('motion_notify_event', self.motion_notify_callback)
        self.canvas = canvas

    def _setEdge(self, **kwargs):
        self._dec_range = np.linspace(-90, 90, self._resolution)
        self._ra_range = np.linspace(-180, 180, self._resolution) + self.proj.ra_0

        # styling: frame needs to be on top of everything, must be transparent
        facecolor = 'None'
        zorder = 1000
        lw = kwargs.pop('lw', 0.7)
        edgecolor = kwargs.pop('edgecolor', 'k')
        # if there is facecolor: clone the polygon and put it in as bottom layer
        facecolor_ = kwargs.pop('facecolor', '#dddddd')

        # polygon of the map edge: top, left, bottom, right
        # don't draw poles if that's a single point
        lines = [self._getMeridian(self.proj.ra_0 + 180, reverse=True), self._getMeridian(self.proj.ra_0 - 180)]
        if not self.proj.poleIsPoint[-90]:
            lines.insert(1, self._getParallel(-90, reverse=True))
        if not self.proj.poleIsPoint[90]:
            lines.insert(0, self._getParallel(90))
        xy = np.concatenate(lines, axis=1).T
        self._edge = Polygon(xy, closed=True, edgecolor=edgecolor, facecolor=facecolor, lw=lw, zorder=zorder,gid="edge", **kwargs)
        self.ax.add_patch(self._edge)

        if facecolor_ is not None:
            zorder = -1000
            edgecolor = 'None'
            poly = Polygon(xy, closed=True, edgecolor=edgecolor, facecolor=facecolor_, zorder=zorder, gid="edge-background")
            self.ax.add_patch(poly)

    def __init__(self, ax, Inpimg, Mask, max_ds=10,PolyAtStart = [(50,50),(100,50),(100,100),(50,100)]):
        self.showverts = True
        self.max_ds = max_ds
        self.Mask = Mask
        self.img = Inpimg
        self.maskinvert = False
        # imshow the image
        self.imgplot = ax.imshow(np.ma.filled(np.ma.array(self.img,mask=self.Mask,fill_value=np.nan)), cmap=cm.gray)
         
        self.poly = Polygon(PolyAtStart, animated=True,
                            fc='y', ec='none', alpha=0.5)
 
        ax.add_patch(self.poly)
        ax.set_clip_on(False)
        ax.set_title("Click and drag a point to move it; "
                     "'i' to insert; 'd' to delete.\n"
                     "'n' to invert the region for masking, 'c' to confirm & apply the mask.")
        self.ax = ax
         
        x, y = zip(*self.poly.xy)
        self.line = Line2D(x, y, color='none', marker='o', mfc='r',
                               alpha=0.7, animated=True)
#        self._update_line()
        self.ax.add_line(self.line)
         
        self.poly.add_callback(self.poly_changed)
        self._ind = None # the active vert
         
        canvas = self.poly.figure.canvas
        canvas.mpl_connect('draw_event', self.draw_callback)
        canvas.mpl_connect('button_press_event', self.button_press_callback)
        canvas.mpl_connect('button_release_event', self.button_release_callback)
        canvas.mpl_connect('key_press_event', self.key_press_callback)
        canvas.mpl_connect('motion_notify_event', self.motion_notify_callback) 
        self.canvas = canvas

    def __init__(self, axtmp, pntxy):
        QtCore.QObject.__init__(self)
        self.ax = axtmp
        self.poly = Polygon([(1, 1)], animated=True)
        self.ax.add_patch(self.poly)
        self.canvas = self.poly.figure.canvas
        self.poly.set_alpha(0.5)
        self.pntxy = pntxy
        self.ishist = True
        self.background = self.canvas.copy_from_bbox(self.ax.bbox)

        xtmp, ytmp = list(zip(*self.poly.xy))

        self.line = Line2D(xtmp, ytmp, marker='o', markerfacecolor='r',
                           color='y', animated=True)
        self.ax.add_line(self.line)

        self.poly.add_callback(self.poly_changed)
        self._ind = None  # the active vert

        self.canvas.mpl_connect('button_press_event',
                                self.button_press_callback)
        self.canvas.mpl_connect('button_release_event',
                                self.button_release_callback)
        self.canvas.mpl_connect('motion_notify_event',
                                self.motion_notify_callback)

 def _setup_patch(self):
     self._patch = Polygon(np.array(list(zip([0, 1], [0, 1]))))
     self._patch.set_zorder(100)
     self._patch.set(**self.plot_opts)
     self._axes.add_patch(self._patch)
     self._patch.set_visible(False)
     self._sync_patch()

def getCatalog(size=10000, surveyname=None):
    # dummy catalog: uniform on sphere
    # Marsaglia (1972)
    xyz = np.random.normal(size=(size, 3))
    r = np.sqrt((xyz**2).sum(axis=1))
    dec = np.arccos(xyz[:,2]/r) / skm.DEG2RAD - 90
    ra = - np.arctan2(xyz[:,0], xyz[:,1]) / skm.DEG2RAD

    if surveyname is not None:
        from matplotlib.patches import Polygon
        # construct survey polygon
        ra_fp, dec_fp = skm.survey_register[surveyname].getFootprint()
        poly = Polygon(np.dstack((ra_fp,dec_fp))[0], closed=True)
        inside = [poly.get_path().contains_point(Point(ra_,dec_)) for (ra_,dec_) in zip(ra,dec)]
        ra = ra[inside]
        dec = dec[inside]

    return ra, dec

    def pressContour(self, event):
        if event.inaxes==None: return
        if event.button != 1: return
        
        # new contour
        if self.poly is None:
            self.poly = Polygon( [(event.xdata , event.ydata)] , animated=False , alpha = .3 , color = 'g')
            self.ax.add_patch(self.poly)
            self.line, = self.ax.plot([event.xdata] , [event.ydata] , 
                                    color = 'g',
                                    linewidth = 2 ,
                                    marker = 'o' ,
                                    markerfacecolor='g', 
                                    animated=False)
            #~ self.canvas.draw()
            self.redraw()
            return
            
        

        # event near a point
        xy = asarray(self.poly.xy)
        xyt = self.poly.get_transform().transform(xy)
        xt, yt = xyt[:, 0], xyt[:, 1]
        
        print '#######'
        d = sqrt((xt-event.x)**2 + (yt-event.y)**2)
        indseq = nonzero(equal(d, amin(d)))[0]
        self._ind = indseq[0]
        if d[self._ind]>=self.epsilon:
            self._ind = None

        self.a=list(numpy.copy(self.poly.xy))
        for i,j in enumerate(self.a):
            self.a[i]=tuple(j)
        
        # new point
        if self._ind is None:
            b=float(event.xdata)
            c=float(event.ydata)
            d=(b,c)
            e=[d]
            #self.a=numpy.append(self.a,e)
            self.a.extend(e)
            self.line.set_xdata( list(self.line.get_xdata()) + [ event.xdata] )
            self.line.set_ydata( list(self.line.get_ydata()) + [ event.ydata] )
            
            #~ self.canvas.draw()
            self.redraw()

        #self.poly.xy=a
        print self.a, type(self.a)
        if self.a != None:    
            test=Path(self.a)
            self.actualSelection, = where(test.contains_points(dot( self.data, self.projection )))
            self.emit(SIGNAL('selectionChanged'))

 def __init__(self, ax):
     """
     :param ax: A matplotlib Axes object to attach the graphical ROI to
     """
     AbstractMplRoi.__init__(self, ax)
     self.plot_opts = {'edgecolor': PATCH_COLOR, 'facecolor': PATCH_COLOR,
                       'alpha': 0.3}
     self._patch = Polygon(np.array(zip([0, 1], [0, 1])))
     self._patch.set_zorder(100)
     self._patch.set(**self.plot_opts)
     self._setup_patch()

    def __init__(self, axes, roi=None):

        super(MplPolygonalROI, self).__init__(axes, roi=roi)

        self.plot_opts = {'edgecolor': PATCH_COLOR,
                          'facecolor': PATCH_COLOR,
                          'alpha': 0.3}

        self._patch = Polygon(np.array(list(zip([0, 1], [0, 1]))), zorder=100)
        self._patch.set_visible(False)
        self._axes.add_patch(self._patch)

    def draw_new_poly(self):
        coords = default_vertices(self.ax)
        Poly = Polygon(coords, animated=True,
                            fc="white", ec='none', alpha=0.2, picker=True)
        self.polyList.append(Poly)
        self.ax.add_patch(Poly)
        x, y = zip(*Poly.xy)
        #line_color = 'none'
        color = np.array((1,1,1))
        marker_face_color = (1,1,1)
        line_width = 4;
        
        line_kwargs = {'lw': line_width, 'color': color, 'mfc': marker_face_color}
        self.line.append(plt.Line2D(x, y, marker='o', alpha=1, animated=True, **line_kwargs))
        self._update_line()
        self.ax.add_line(self.line[-1])

        Poly.add_callback(self.poly_changed)
        Poly.num = self.poly_count
        self._ind = None  # the active vert
        self.poly_count+=1;

    def __init__(self, filename, **kwargs):
        """
        Create a new coastline polygon.

        Parameters
        ----------
        color : color, optional, default 'gray'
            Line color of the coastline.
        land : color, optional, default 'seashell'
            Fill color of the land polygons.

        Other parameters
        ----------------
        kwargs : polygon properties
            Other parameters passed on to :class:`~matplotlib.patches.Polygon`,
            e.g. zorder=N to control drawing the land polygons above/below
            other data.
        """
        color = kwargs.pop('color', 'gray')
        land = kwargs.pop('land', 'seashell')
        self.data = []
        self.extents = None

        if not color:
            color = 'none'
        if not land:
            land = 'none'

        xy = [[None, None], [None, None]]
        Polygon.__init__(self, xy, edgecolor=color, facecolor=land, **kwargs)

        datapath = os.path.join(os.path.dirname(__file__), 'data')
        coastfile = _find(filename, datapath)
        if coastfile:
            file = shapefile.Reader(coastfile)
            for shape in file.shapes():
                for points in _split(shape.points, shape.parts):
                    self.data += [Path(points)]
        else:
            raise Warning('coastline "%s" not found in directory "%s"' % (filename, datapath))