def number_of_images(self, sourceposition):
        rc=self.radial_caustic()
        tc=self.tangential_caustic()

        if usePath:
        
            # New Matplotlib:
            rc=Path(rc)
            tc=Path(tc)
            if rc.contains_points(np.atleast_2d(sourceposition))[0]:
               if tc.contains_points(np.atleast_2d(sourceposition))[0]:
                   return 4
               return 2
            if tc.contains_points(np.atleast_2d(sourceposition))[0]:
               return 3
        
        else:
            # Old Matplotlib:
            if nxutils.points_inside_poly(np.atleast_2d(sourceposition),rc)[0]:
                if nxutils.points_inside_poly(np.atleast_2d(sourceposition),tc)[0]:
                    return 4
                return 2
            if nxutils.points_inside_poly(np.atleast_2d(sourceposition),tc)[0]:
                return 3

        return 1

 def number_of_images(self, sourceposition):
     if nx.points_inside_poly(np.atleast_2d(sourceposition), self.radial_caustic())[0]:
         if nx.points_inside_poly(np.atleast_2d(sourceposition), self.tangential_caustic())[0]:
             return 4
         return 2
     if nx.points_inside_poly(np.atleast_2d(sourceposition), self.tangential_caustic())[0]:
         return 3
     return 1

def main(shapefile, picklefile):
    if picklefile:
      [npts, x, y, z, zraw, xil, yil, grid, missing] = cPickle.load(open(picklefile,'rb'))
      points = np.vstack((x,y)).T

    # Load administrative area
    shp = ShapeFile(shapefile)
    dbf = dbflib.open(shapefile)

    coastline = []

    # Process every shape from the ShapeFile
    print "Processing shapes ..."
    for npoly in range(shp.info()[0]):
        shp_object = shp.read_object(npoly)
        shp_dict = dbf.read_record(npoly)
        verts = shp_object.vertices()

        if "NAME_1" in shp_dict:
          name = "%s" % (shp_dict["NAME_1"])
        else:
          name = "Unknown"

        print "Processing %s" % (name)
        # Extract city polygon vertices (ring per ring)
        for ring in verts:
          vx = []
          vy = []
          for point in ring:
            vx.append(point[0])
            vy.append(point[1])

          # Only process big enough rings
          if len(vx) > 256: # big enough
            poly_verts = zip(vx,vy)

            if picklefile:
              # Compute intersections with the city
              intersection = points_inside_poly(points, poly_verts)
              npts = sum(1 for x in points_inside_poly(points, poly_verts) if x)
            else:
              npts = 1 # Add this polygon

            # Add the ring to the coastine if measures inside
            if npts > 0:
              polygon = Polygon(poly_verts)
              if not polygon.is_empty and polygon.is_valid:
                print "- Add polygon (%d)" % (len(vx))
                coastline.append(polygon)
            else:
                print "- Skip polygon (%d)" % (len(vx))
    
    print "Union of %d polygons" % len(coastline)
    coast = cascaded_union(coastline)
    cPickle.dump(coast,open('coastline.pickle','wb'),-1)
    print "Done."

def get_contour_mask(doselut, dosegridpoints, contour):
    """Get the mask for the contour with respect to the dose plane."""

    grid = nx.points_inside_poly(dosegridpoints, contour)
    grid = grid.reshape((len(doselut[1]), len(doselut[0])))

    return grid

def inregion(x,y,r):
    """
    ins = inregion(x,y,r)

    Returns an array of booleans indicating whether the x,y pairs are
    in region r. If region r contains multiple polygons, the ones inside 
    the biggest one are assumed to be holes; the ones inside holes
    are assumed to be islands; and so on.
    
    :Parameters:
      x : array
        x coordinates of test points
      y : array
        y coordinates of test points
      r : ShapeObject
        The region.
    """
    xy = np.vstack((x,y)).T

    # Record whether each point is inside each polygon.
    ins = []
    for v in r:
        ins.append(points_inside_poly(xy,v))
    ins = np.array(ins)
    
    # Return an array of booleans. An element is True if
    # the corresponding point is inside an odd number of polygons.
    return np.sum(ins, axis=0) % 2 == 1

def points_in_polys(points, polys):
	# Function to quickly check stranding of many points
	insidePoly = np.array([False]*len(points))
	for poly in polys:
		# NB: use contains_points for matplotlib version >= 1.2.0
		insidePoly[nx.points_inside_poly(points, poly)] = True
	return insidePoly

def inside_poly(vertices,data):
    if(matplotlib.__version__ < '1.2'):
        mask = points_inside_poly(data, vertices)
    else:
        mask = Path(vertices).contains_points(data)

    return mask

	def __find(self, verts):
		ps = np.array(map(lambda x: [x[1], x[2]], self.points), float)
		points_in = nx.points_inside_poly(ps, verts)
		for i in range(len(points_in)):
			if points_in[i] == True:
				 return self.points[i]
		return None

def poly_over_under(ra,dec,patch_left,patch_right):
    '''
    Short cut to Polygon.Polygon
    (p and q are polygons)
    p & q: intersection: a polygon with the area that is covered by both p and q
    p | q: union: a polygon containing the area that is covered by p or q or both
    p - q: difference: a polygon with the area of p that is not covered by q
    p + q: sum: same as union
    p ^ q: xor: a polygon with the area that is covered by exactly one of p and q

    len(p):
    p[i]:

    number of contours
    contour with index i, the same as p.contour(i), slicing is not yet supported
    '''
    # returns indices of ra, dec that are not in the overlap region 
    # of the two patches.
    # use ra,dec.
    lr,ld = mpl_patch_XY2radec(patch_left)
    rr,rd = mpl_patch_XY2radec(patch_right)
    # Polygon needs tuples...
    left = Polygon.Polygon(ndarray2tuple(np.column_stack((lr,ld))))
    right = Polygon.Polygon(ndarray2tuple(np.column_stack((rr,rd))))
    # p & q: intersection: a polygon with the area that is covered by both p and q
    overlap = right & left
    # vertices of overlapped region
    if len(overlap[0]) > 0:
        verts = np.transpose(np.array(overlap[0]))
        radec = np.column_stack((ra,dec))
        mask = nxutils.points_inside_poly(radec, np.array(overlap[0]))
    # sloppy? flip T/F: subtract 1 from bool array and then get rid of neg sign
        not_overlapped = np.nonzero(abs(mask-1))[0]
    else: not_overlapped = []
    return not_overlapped

def _raster_points(tmpx, tmpy, gridShape):
    """
    Find the raster grid points that lie within the voxel
    """
    if max(tmpx) < 0 or max(tmpy) < 0 or min(tmpx) >= gridShape[1] or min(tmpy) >= gridShape[0]:
        # points lie outside the rasterization grid
        # so, none of them are good.
        return ([], [])

    resVol = zip(tmpx[[0, 1, 2, 3, 0]], tmpy[[0, 1, 2, 3, 0]])

    # Getting all of the points that the polygon has, and then some.
    # This meshed grid is bounded by the domain.
    bbox = (
        (int(max(np.floor(min(tmpy)), 0)), int(min(np.ceil(max(tmpy)), gridShape[0] - 1))),
        (int(max(np.floor(min(tmpx)), 0)), int(min(np.ceil(max(tmpx)), gridShape[1] - 1))),
    )
    (ygrid, xgrid) = np.meshgrid(np.arange(bbox[0][0], bbox[0][1] + 1), np.arange(bbox[1][0], bbox[1][1] + 1))
    gridPoints = zip(xgrid.flat, ygrid.flat)

    if len(gridPoints) == 0:
        print("Bad situation...:", bbox, gridShape, min(tmpy), max(tmpy), min(tmpx), max(tmpx))
        gridPoints = np.zeros((0, 2), dtype="i")

    # Determines which points fall within the resolution volume.  These
    # points will be the ones that will be assigned the value of the
    # original data point that the resolution volume represents.
    goodPoints = points_inside_poly(gridPoints, resVol)

    return (ygrid.flat[goodPoints], xgrid.flat[goodPoints])

    def getMDAreas(self,dx=20000,dy=20000,
                    llcrnrlon=-119.2,
                    llcrnrlat=23.15,
                    urcrnrlon=-65.68,
                    urcrnrlat=48.7):
        bmap = Basemap(projection="lcc",
                       llcrnrlon=llcrnrlon,
                       llcrnrlat=llcrnrlat,
                       urcrnrlon=urcrnrlon,
                       urcrnrlat=urcrnrlat,
                       resolution='l',
                       lat_0=38.5,
                       lat_1=38.5,
                       lon_0=-97.0)
        from matplotlib.nxutils import points_inside_poly
        xs = np.arange(bmap.llcrnrx,bmap.urcrnrx + dx,dx)
        ys = np.arange(bmap.llcrnry,bmap.urcrnry + dy,dy)
        lon,lat,x_grid,y_grid = bmap.makegrid(xs.shape[0],ys.shape[0],returnxy=True)
        x, y = x_grid.flatten(), y_grid.flatten()
        points = np.vstack((x,y)).T
        
        nx = xs.shape[0]
        ny = ys.shape[0]
        areas = np.zeros((self.data.shape[0],))

        for i in xrange(self.data.shape[0]):
            md_x,md_y = bmap(self.data['Lon'][i],self.data['Lat'][i])
            poly_xy = np.vstack((md_x,md_y)).T
            areas[i] = np.nonzero(points_inside_poly(points,poly_xy))[0].shape[0] * dx * dy / 1000**2

        return areas

 def mask_polygon(self, polyverts, mask_value=0.0):
     """
     Mask Cartesian points contained within the polygon defined by polyverts
     
     A cell is masked if the cell center (x_rho, y_rho) is within the
     polygon. Other sub-masks (mask_u, mask_v, and mask_psi) are updated
     automatically.
     
     mask_value [=0.0] may be specified to alter the value of the mask set
     within the polygon.  E.g., mask_value=1 for water points.
     """
     
     polyverts = np.asarray(polyverts)
     assert polyverts.ndim == 2, \
         'polyverts must be a 2D array, or a similar sequence'
     assert polyverts.shape[1] == 2, \
         'polyverts must be two columns of points'
     assert polyverts.shape[0] > 2, \
         'polyverts must contain at least 3 points'
     
     mask = self.mask_rho
     inside = points_inside_poly(
         np.vstack( (self.x_rho.flat, self.y_rho.flat) ).T,
         polyverts)
     if np.any(inside):
         self.mask_rho.flat[inside] = mask_value

    def get_variables(self, requestedVariables, time=None,
                      x=None, y=None, z=None, block=False):

        if isinstance(requestedVariables, str):
            requestedVariables = [requestedVariables]

        self.check_arguments(requestedVariables, time, x, y, z)

        # Apparently it is necessary to first convert from x,y to lon,lat
        # using proj library, and back to x,y using Basemap instance
        # Perhaps a bug in Basemap related to x_0/y_0 offsets?
        # Nevertheless, seems not to affect performance
        lon, lat = self.xy2lonlat(x, y)
        x, y = self.map(lon, lat, inverse=False)
        points = np.c_[x, y]

        insidePoly = np.array([False]*len(x))
        if has_nxutils is True:
            for polygon in self.polygons:
                insidePoly[nx.points_inside_poly(points, polygon)] = True
        else:
            for polygon in self.polygons:
                insidePoly += np.array(polygon.contains_points(points))

        variables = {}
        variables['land_binary_mask'] = insidePoly

        return variables

    def contains(self, x, y):
        """
        Test whether a set of (x,y) points falls within
        the region of interest

        :param x: A list of x points
        :param y: A list of y points

        *Returns*

           A list of True/False values, for whether each (x,y)
           point falls within the ROI

        """
        if not self.defined():
            raise UndefinedROI
        if not isinstance(x, np.ndarray):
            x = np.asarray(x)
        if not isinstance(y, np.ndarray):
            y = np.asarray(y)

        xypts = np.column_stack((x.flat, y.flat))
        xyvts = np.column_stack((self.vx, self.vy))
        result = points_inside_poly(xypts, xyvts)
        good = np.isfinite(xypts).all(axis=1)
        result[~good] = False
        result.shape = x.shape
        return result

 def calculate_polygon_intersection (self, geometry, points):
     ''' Apply nxutils algorithm to determine point-in-polygon relationship. '''
     rings = geometry ['coordinates']
     ring_points = []
     for ring in rings:
         for point in ring:
             ring_points.append (list (point))
     numpy_geometry = numpy.array (ring_points)
     return nx.points_inside_poly (points, numpy_geometry)

def inpolygon(x,y,xp,yp):
  '''
  Points inside polygon test.

  Based on matplotlib nxutils for old matplotlib versions
  http://matplotlib.sourceforge.net/faq/howto_faq.html#test-whether-a-point-is-inside-a-polygon

  Can also use Path.contains_point, for recent versions, or
  the pnpoly extension - point in polyon test - by W R Franklin,
  http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
  '''

  try:
    from matplotlib.nxutils import pnpoly, points_inside_poly
    _use_mpl = 1
  except ImportError:
    from matplotlib.path import Path
    _use_mpl = 2
  except:
    import pnpoly
    _use_mpl = False

  shape=False
  try:
    if x.ndim>1:
      shape=x.shape
      x=x.flat
      y=y.flat
  except: pass

  if _use_mpl==1:
    verts=np.array(zip(xp,yp))
    if isiterable(x,y):
      points=np.array(zip(x,y))
      res=points_inside_poly(points,verts)
    else:
      res=pnpoly(x,y,verts)==1

  elif _use_mpl==2:
    verts=np.array(zip(xp,yp))
    p=Path(verts)
    if not isiterable(x,y): x,y,itr=[x],[y],0
    else: itr=1
    res=[p.contains_point((x[i],y[i]), radius=0.) for i in range(len(x))]
    res=np.asarray(res,'bool')
    if not itr: res=res[0]

  else:
    if not isiterable(x,y): x,y,itr=[x],[y],0
    else: itr=1
    res=np.zeros(len(x))
    res=[pnpoly.pnpoly(x[i],y[i],xp,yp,len(xp)) for i in range(len(x))]
    res=np.asarray(res)==1
    if not itr: res=res[0]

  if not shape is False: res.shape=shape
  return res

    def outer_loop(G, cycles):
        coms = array([c.com for c in cycles])

        for c in cycles:
            xyverts = zeros((c.coords.shape[0] + 1, c.coords.shape[1]))
            xyverts[:-1, :] = c.coords
            xyverts[-1, :] = c.coords[0, :]
            if nxutils.points_inside_poly(coms, xyverts).all():
                return c

 def getROI(self, verts):
     ind = points_inside_poly(self.xys, verts)
     self.canvas.draw_idle()
     self.canvas.widgetlock.release(self.lasso)
     del self.lasso
     self.ROI = ind
     self.ROI = self.ROI.reshape((600,800), order='F')
     self.canvas.mpl_disconnect(self.cid)
     self.onDraw()

 def filter(self):
     """ Data flow into here """
     while True:
         a = (yield)
         # good = np.ones(a.shape, dtype=bool)
         
         coord0 = self.coord_names[0]
         coord1 = self.coord_names[1]
         in_poly_mask = points_inside_poly(zip(a[coord0], a[coord1]), self.verts) == 1                    
         self.target.send(a[in_poly_mask])

def dpoly(p,pv):
    """Signed distance function for polygon with vertices pv.

    Usually pv should also be provided as fixed points in distmesh2d.

    pv should be provided as a list of coordinates [(x0,y0), (x1,y1), ...]
    or an array of shape (nv, 2).
    """
    from matplotlib.nxutils import points_inside_poly
    return (-1)**points_inside_poly(p, pv) * dsegment(p, pv).min(1)