def doPolygonize():
  blocks = polygonize(lines)
  writeBlocks(blocks, args[0] + '-blocks.geojson')

  blocks = polygonize(lines)
  bounds = Polygon([
    [minlng, minlat],
    [minlng, maxlat],
    [maxlng, maxlat],
    [maxlng, minlat],
    [minlng, minlat]
  ])
  # Geometry transform function based on pyproj.transform
  project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:3785'),
    pyproj.Proj(init='EPSG:4326'))
  print bounds
  print transform(project, bounds)

  print 'finding holes'
  for index, block in enumerate(blocks):
    if index % 1000 == 0:
      print "diff'd  %s" % (index)
    if not block.is_valid:
      print explain_validity(block)
      print transform(project, block)
    else:
      bounds = bounds.difference(block)
  print bounds

def create_union(in_ply1, in_ply2, result_geojson):
    """
    Create union polygon
    :param in_ply1: first input shapely polygon
    :param in_ply2: second input shapely polygon
    :param result_geojson: output geojson file including full file path
    :return: shapely MultiPolygon
    """
    # union the polygon outer linestrings together
    outer_bndry = in_ply1.boundary.union(in_ply2.boundary)

    # rebuild linestrings into polygons
    output_poly_list = polygonize(outer_bndry)

    out_geojson = dict(type='FeatureCollection', features=[])

    # generate geojson file output
    for (index_num, ply) in enumerate(output_poly_list):
        feature = dict(type='Feature', properties=dict(id=index_num))
        feature['geometry'] = ply.__geo_interface__
        out_geojson['features'].append(feature)

    # create geojson file on disk
    json.dump(out_geojson, open(result_geojson, 'w'))

    # create shapely MultiPolygon
    ply_list = []
    for fp in polygonize(outer_bndry):
        ply_list.append(fp)

    out_multi_ply = MultiPolygon(ply_list)

    return out_multi_ply

def alpha_shape(points, alpha):
    def add_edge(points, i, j):
        if (i, j) in edges or (j, i) in edges:
            return
        edges.add((i, j))
        edge_points.append((points[i], points[j]))
    tri = DelaunayTri(points)
    edges = set()
    edge_points = []
    for i1, i2, i3 in tri.vertices:
        x1, y1 = tri.points[i1]
        x2, y2 = tri.points[i2]
        x3, y3 = tri.points[i3]
        a = hypot(x1 - x2, y1 - y2)
        b = hypot(x2 - x3, y2 - y3)
        c = hypot(x3 - x1, y3 - y1)
        s = (a + b + c) / 2.0
        area = sqrt(s * (s - a) * (s - b) * (s - c))
        radius = a * b * c / (4 * area)
        if radius < 1.0 / alpha:
            add_edge(tri.points, i1, i2)
            add_edge(tri.points, i2, i3)
            add_edge(tri.points, i3, i1)
    shape = cascaded_union(list(polygonize(MultiLineString(edge_points))))
    return shape

def alpha_shape(points, alpha):
    from shapely.ops import cascaded_union, polygonize
    from scipy.spatial import Delaunay
    import numpy as np
    import math
    import shapely.geometry as geometry
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    """
    if len(points) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull

    def add_edge(edges, edge_points, coords, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
                # already added
            return
        edges.add( (i, j) )
        edge_points.append(coords[ [i, j] ])
    coords = np.array([point.coords[0]
                       for point in points])
    tri = Delaunay(coords)
    edges = set()
    edge_points = []
    # loop over triangles:
    # ia, ib, ic = indices of corner points of the
    # triangle
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]
        # Lengths of sides of triangle
        a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
        b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
        c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
        # Semiperimeter of triangle
        s = (a + b + c)/2.0
        # Area of triangle by Heron's formula
        area = math.sqrt(s*(s-a)*(s-b)*(s-c))
        circum_r = a*b*c/(4.0*area)
        # Here's the radius filter.
        #print circum_r
        if circum_r < 1.0/alpha:
            add_edge(edges, edge_points, coords, ia, ib)
            add_edge(edges, edge_points, coords, ib, ic)
            add_edge(edges, edge_points, coords, ic, ia)
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return cascaded_union(triangles), edge_points

    def get_boundingpolygon(self):
        """
            TODO: Implement ncell bbox

            order of lines creation (assumes 0,0 is x)
            -----3-----
            |         |
            4         2
            |         |
            x----1-----
        """

        if self._ndim == 2: # CGRID
            nx,ny = self._xarray.shape
            one = MultiLineString([((self._xarray[i][0],self._yarray[i][0]),(self._xarray[i+1][0],self._yarray[i+1][0])) for i in range(nx-1)])
            two = MultiLineString([((self._xarray[nx-1][j],self._yarray[nx-1][j]),(self._xarray[nx-1][j+1],self._yarray[nx-1][j+1])) for j in range(ny-1)])
            three = MultiLineString([((self._xarray[i][ny-1],self._yarray[i][ny-1]),(self._xarray[i-1][ny-1],self._yarray[i-1][ny-1])) for i in reversed(list(range(1,nx)))])
            four = MultiLineString([((self._xarray[0][j],self._yarray[0][j]),(self._xarray[0][j-1],self._yarray[0][j-1])) for j in reversed(list(range(1,ny)))])
            m = one.union(two).union(three).union(four)
        else: # RGRID
            nx,ny = self._xarray.shape[0], self._yarray.shape[0]
            one = LineString([(self._xarray[i], self._yarray[0]) for i in range(nx)])
            two = LineString([(self._xarray[-1], self._yarray[i]) for i in range(ny)])
            three = LineString([(self._xarray[i], self._yarray[-1]) for i in reversed(list(range(nx)))])
            four = LineString([(self._xarray[0], self._yarray[i]) for i in reversed(list(range(ny)))])
            m = MultiLineString([one,two,three,four])

        polygons = list(polygonize(m))
        # -- polygonize returns a list of polygons, including interior features, the largest in area "should" be the full feature
        assert len(polygons) > 0, "Could not determine a polygon"
        polygon = sorted(polygons, key=lambda x: x.area)[-1]
        return polygon

def alpha_shape(x, y, alpha):
    coords = np.c_[x, y]
    tri = scipy.spatial.Delaunay(coords)
    
    # edges = []
    edge_coords = []
    for ia, ib, ic in tri.vertices:
        ab = math.sqrt((coords[ia,0]-coords[ib,0])**2 + (coords[ia,1]-coords[ib,1])**2)
        bc = math.sqrt((coords[ib,0]-coords[ic,0])**2 + (coords[ib,1]-coords[ic,1])**2)
        ac = math.sqrt((coords[ia,0]-coords[ic,0])**2 + (coords[ia,1]-coords[ic,1])**2)
        
        semiperim = 0.5*(ab+bc+ac)
        area = math.sqrt(semiperim*(semiperim-ab)*(semiperim-bc)*(semiperim-ac))
        radius = 0.25*ab*bc*ac/area
        
        if radius < 1.0/alpha:
            # edges.append((tuple(coords[ia]), tuple(coords[ib])))
            # edges.append((tuple(coords[ib]), tuple(coords[ic])))
            # edges.append((tuple(coords[ia]), tuple(coords[ic])))
            edge_coords.append(coords[[ia, ib]])
            edge_coords.append(coords[[ib, ic]])
            edge_coords.append(coords[[ia, ic]])
            
    # edges = set(edges)
    m = MultiLineString(edge_coords)
    return cascaded_union(list(polygonize(m))), edge_coords

def difference(line1, line2, close_ends=False):
    """ Create polygons from two LineString objects.

    Parameters:
        line1 (LineString) : a line representing the initial condition.
        line2 (LineString) : a line representing the final condition.
        close_ends (bool) : option to close open line ends with vertical line segments.

    Returns:
        intersections (Point array) : the intersections between the LineString objects.
        polygons (Polygon array) : the polygons between the lines.
        signs (int array) : contains values of +1 or -1 to identify polygons as cut or fill.

    """
    if close_ends==True:
        line1, line2 = close(line1, line2)

    intersections = line1.intersection(line2)

    segs1 = cut_by_points([line1], intersections)
    segs2 = cut_by_points([line2], intersections)

    polygons = polygonize([segs1, segs2])

    signs = sign(linemerge(segs1), linemerge(segs2))

    # can't pass the polygonize generator to my class so convert the polygons into an array
    polygontxt = []
    areas = []
    for i, poly in enumerate(polygons):
        polygontxt.append(poly)
        areas.append(poly.area*signs[i])
    cutfill = pnd.Series(asarray(areas), name='area')

    return intersections, polygontxt, cutfill

    def _cont_to_polys(self, cont_lats, cont_lons):
        polys = []

        start_idx = 0
        splits = []
        while True:
            try:
                split_idx = cont_lats[start_idx:].index(99.99)
                splits.append(start_idx + split_idx)
                start_idx += split_idx + 1
            except ValueError:
                break

        splits = [ -1 ] + splits + [ len(cont_lats) + 1 ]
        poly_lats = [ cont_lats[splits[i] + 1:splits[i + 1]] for i in xrange(len(splits) - 1) ]
        poly_lons = [ cont_lons[splits[i] + 1:splits[i + 1]] for i in xrange(len(splits) - 1) ]

        # Intersect with the US boundary shape file.
        for plat, plon in zip(poly_lats, poly_lons):
            cont = LineString(zip(plon, plat))

            if plat[0] != plat[-1] or plon[0] != plon[-1]:
                # If the line is not a closed contour, then it intersects with the edge of the US. Extend
                #   the ends a little bit to make sure it's outside the edge.
                dln = np.diff(plon)
                dlt = np.diff(plat)

                pre = [ (plon[0] - 0.5 * dln[0], plat[0] - 0.5 * dlt[0]) ]
                post = [ (plon[-1] + 0.5 * dln[-1], plat[-1] + 0.5 * dlt[-1]) ]

                cont.coords = pre + list(cont.coords) + post 

            # polygonize() will split the country into two parts: one inside the outlook and one outside.
            #   Construct test_ln that is to the right of (inside) the contour and keep only the polygon
            #   that contains the line
            test_ln = cont.parallel_offset(0.05, 'right')

            for poly in polygonize(self._conus.boundary.union(cont)):
                if (poly.crosses(test_ln) or poly.contains(test_ln)) and self._conus.contains(poly.buffer(-0.01)):
                    polys.append(poly)

        # Sort the polygons by area so we intersect the big ones with the big ones first.
        polys.sort(key=lambda p: p.area, reverse=True)

        # If any polygons intersect, replace them with their intersection.
        intsct_polys = []
        while len(polys) > 0:
            intsct_poly = polys.pop(0)
            pops = []
            for idx, poly in enumerate(polys):
                if intsct_poly.intersects(poly):
                    intsct_poly = intsct_poly.intersection(poly)
                    pops.append(idx)

            for pop_idx in pops[::-1]:
                polys.pop(pop_idx)

            intsct_polys.append(intsct_poly)
        return intsct_polys

def procesaLineaInterna(featuresExternas, featuresInternas, featuresCentroide, featureDefn):
	#print "Procedemos a procesar las lineas internas"
	
	centroides = []
	for centroide in featuresCentroide:
		#obtenemos la altura y el rotulo del estilo de cada centroide
		for n in centroide.GetStyleString().split(','):
			if n.startswith('s'):
				altura = float(n.replace('s:', '').replace('g', ''))
			elif n.startswith('t'):
				rotulo = n.split('"')[1]
		punto = centroide.GetGeometryRef()
		x = punto.GetX()
		y = punto.GetY()
		longitudRotulo = len(rotulo)
		factor = 0.15 * (altura * 3.3333)
		desfaseX = longitudRotulo * factor - 0.05
		punto.SetPoint(point = 0, x = x + desfaseX, y = y - 0.20)
		
		centroides.append((rotulo, punto))
	
	featuresProceso = featuresExternas + featuresInternas
	
	outFeature = []
	if len(featuresProceso) > 1:
		geometry_out = None
		for inFeature in featuresProceso:
			geometry_in = inFeature.GetGeometryRef()
			if geometry_out is None:
				geometry_out = geometry_in
				geometry_out = ogr.ForceToMultiLineString(geometry_out)
			else:
				geometry_out = geometry_out.Union(geometry_in) 
		
		lineasInternasShapely = loads(geometry_out.ExportToWkt())
		polygonsShapely = polygonize(lineasInternasShapely)
	
		polygonGeom = []
		for polygon in polygonsShapely:
			polygonGeom.append(ogr.CreateGeometryFromWkt(dumps(polygon)))
		
		for pol in polygonGeom:
			for cen in centroides:
				if pol.Contains(cen[1]):
					feature = ogr.Feature(featureDefn)
					feature.SetGeometry(pol)
					feature.SetField('rotulo', cen[0])
					outFeature.append(feature.Clone())
					feature.Destroy()
	else:
		feature = ogr.Feature(featureDefn)
		geometryPoly = ogr.BuildPolygonFromEdges(ogr.ForceToMultiLineString(featuresProceso[0].GetGeometryRef()), dfTolerance = 0)
		feature.SetGeometry(geometryPoly)
		feature.SetField('rotulo', centroides[0][0])
		outFeature.append(feature.Clone())
		feature.Destroy()
	
	
	return outFeature

def get_convex_hull(points):
    """ Find the convex hull of points, return as Shapely polygon """
    log.info("Finding convex hull of %i points", len(points))
    hull = ConvexHull(points)
    log.info("Found convex hull with %i facets", len(hull.simplices))
    hull_edges = []
    for simplex in hull.simplices:
        hull_edges.append(
            zip(points[simplex, 0], points[simplex, 1]))
    polygon = list(polygonize(hull_edges))
    assert(len(polygon) == 1)
    return polygon[0]

def get_shapes(filelike):
    shapes = []
    node_cache = {}
    way_cache = {}
    for thing in iter_osm_file(filelike):
        if type(thing) == Node:
            pt = (thing.lon, thing.lat)
            shape = Point(pt)

            node_cache[thing.id] = pt

            if thing.tags:
                shapes.append((thing, shape))

        elif type(thing) == Way:
            points = []
            for nd in thing.nds:
                node_loc = node_cache.get(nd)
                if node_loc:
                    points.append(node_loc)
                else:
                    raise Exception("Way %s references node %s which is not parsed yet." % (thing.id, nd))

            if way_is_polygon(thing):
                shape = Polygon(points)
            else:
                shape = LineString(points)

            way_cache[thing.id] = points

            if any(thing.tags):
                # Only include tagged things at this point. Otherwise,
                # the shapes that are part of multipolygon relations
                # will be included twice.
                shapes.append((thing, shape))

        elif type(thing) == Relation:
            if any([t.key == 'type' and t.value == 'multipolygon' for t in thing.tags]):
                parts = []
                for member in thing.members:
                    if member.type == 'way':
                        shape = way_cache.get(member.ref)
                        if not shape:
                            raise Exception("Relation %s references way %s which is not parsed yet." % (thing.id, member.ref))

                        parts.append(shape)

                # Polygonize will return all the polygons created, so the
                # inner parts of the multipolygons will be returned twice
                # we only want the first one
                shapes.append((thing, next(polygonize(parts))))

    return shapes

 def test_polygonize(self):
     lines = [
         LineString(((0, 0), (1, 1))),
         LineString(((0, 0), (0, 1))),
         LineString(((0, 1), (1, 1))),
         LineString(((1, 1), (1, 0))),
         LineString(((1, 0), (0, 0))),
         LineString(((5, 5), (6, 6))),
         Point(0, 0),
         ]
     result = list(polygonize(lines))
     self.assertTrue(all([isinstance(x, Polygon) for x in result]))

def save(datasource, filename):
    """ Save a Datasource instance to a named OGR datasource.
    """
    ext = splitext(filename)[1]
    
    out_driver = ogr.GetDriverByName(drivers.get(ext))
    out_source = out_driver.CreateDataSource(filename)
    
    if out_source is None:
        raise Exception('Failed creation of %s - is there one already?' % filename)
    
    out_layer = out_source.CreateLayer('default', datasource.srs, ogr.wkbMultiPolygon)
    
    for field in datasource.fields:
        field_defn = ogr.FieldDefn(field.name, field.type)
        field_defn.SetWidth(field.width)
        out_layer.CreateField(field_defn)
    
    for i in datasource._indexes():
    
        segments = datasource.db.execute("""SELECT x1, y1, x2, y2
                                            FROM segments
                                            WHERE (src1_id = ? OR src2_id = ?)
                                              AND removed = 0""", (i, i))
    
        lines = [datasource.memo_line(x1, y1, x2, y2) for (x1, y1, x2, y2) in segments]
        
        try:
            poly = polygonize(lines).next()
    
        except StopIteration:
            lost_area = datasource.shapes[i].area
            lost_portion = lost_area / (datasource.tolerance ** 2)
            
            if lost_portion < 4:
                # It's just small.
                print >> stderr, 'Skipped small feature #%(i)d' % locals()
                continue
    
            # This is a bug we don't understand yet.
            raise Exception('Failed to get a meaningful polygon out of large feature #%(i)d' % locals())
        
        feat = ogr.Feature(out_layer.GetLayerDefn())
        
        for (j, field) in enumerate(datasource.fields):
            feat.SetField(field.name, datasource.values[i][j])
        
        geom = ogr.CreateGeometryFromWkb(dumps(poly))
        
        feat.SetGeometry(geom)
    
        out_layer.CreateFeature(feat)

def intersect(a, b):
    ap, bp = [orient_offset_polygon(p) for p in [a,b]]
    overlap = ap.intersection(bp)

    if overlap.is_empty:
        return []

    if not hasattr(overlap, 'exterior'):
        return [Poly(p.exterior.coords, p.centroid.coords[0], 0)
                for p in polygonize(overlap)]
    
    return [Poly(overlap.exterior.coords,
                 overlap.centroid.coords[0], 0)]

 def processAlgorithm(self, progress):
     try:
         from shapely.ops import polygonize
         from shapely.geometry import Point, MultiLineString
     except ImportError:
         raise GeoAlgorithmExecutionException(
                 'Polygonize algorithm requires shapely module!')
     vlayer = dataobjects.getObjectFromUri(
             self.getParameterValue(self.INPUT))
     output = self.getOutputFromName(self.OUTPUT)
     vprovider = vlayer.dataProvider()
     if self.getParameterValue(self.FIELDS):
         fields = vprovider.fields()
     else:
         fields = QgsFields()
     if self.getParameterValue(self.GEOMETRY):
         fieldsCount = fields.count()
         fields.append(QgsField('area', QVariant.Double, 'double', 16, 2))
         fields.append(QgsField('perimeter', QVariant.Double,
                                'double', 16, 2))
     allLinesList = []
     features = vector.features(vlayer)
     current = 0
     total = 40.0 / float(len(features))
     for inFeat in features:
         inGeom = inFeat.geometry()
         if inGeom.isMultipart():
             allLinesList.extend(inGeom.asMultiPolyline())
         else:
             allLinesList.append(inGeom.asPolyline())
         current += 1
         progress.setPercentage(int(current * total))
     progress.setPercentage(40)
     allLines = MultiLineString(allLinesList)
     allLines = allLines.union(Point(0, 0))
     progress.setPercentage(45)
     polygons = list(polygonize([allLines]))
     progress.setPercentage(50)
     writer = output.getVectorWriter(fields, QGis.WKBPolygon, vlayer.crs())
     outFeat = QgsFeature()
     current = 0
     total = 50.0 / float(len(polygons))
     for polygon in polygons:
         outFeat.setGeometry(QgsGeometry.fromWkt(polygon.wkt))
         if self.getParameterValue(self.GEOMETRY):
             outFeat.setAttributes([None] * fieldsCount + [polygon.area,
                                   polygon.length])
         writer.addFeature(outFeat)
         current += 1
         progress.setPercentage(50 + int(current * total))
     del writer

def create_union(in_ply1, in_ply2):
    """
    Create union polygon
    :param in_ply1: first input polygon
    :param in_ply2: second input polygon
    :return: shapely polgon
    """
    # union the polygon outer linestrings together
    outer_bndry = in_ply1.boundary.union(in_ply2.boundary)

    # rebuild linestrings into polygons
    output_poly = polygonize(outer_bndry)
    # pprint(list(polygonize(outer_bndry)))
    return output_poly

    def processAlgorithm(self, progress):
        vlayer = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT))
        output = self.getOutputFromName(self.OUTPUT)
        vprovider = vlayer.dataProvider()
        if self.getParameterValue(self.FIELDS):
            fields = vprovider.fields()
        else:
            fields = QgsFields()
        if self.getParameterValue(self.GEOMETRY):
            fieldsCount = fields.count()
            fields.append(QgsField('area', QVariant.Double, 'double', 16, 2))
            fields.append(QgsField('perimeter', QVariant.Double,
                                   'double', 16, 2))
        allLinesList = []
        features = vector.features(vlayer)
        progress.setInfo(self.tr('Processing lines...'))
        total = 40.0 / len(features)
        for current, inFeat in enumerate(features):
            inGeom = inFeat.geometry()
            if inGeom.isMultipart():
                allLinesList.extend(inGeom.asMultiPolyline())
            else:
                allLinesList.append(inGeom.asPolyline())
            progress.setPercentage(int(current * total))

        progress.setPercentage(40)
        allLines = MultiLineString(allLinesList)

        progress.setInfo(self.tr('Noding lines...'))
        allLines = unary_union(allLines)

        progress.setPercentage(45)
        progress.setInfo(self.tr('Polygonizing...'))
        polygons = list(polygonize([allLines]))
        if not polygons:
            raise GeoAlgorithmExecutionException(self.tr('No polygons were created!'))
        progress.setPercentage(50)

        progress.setInfo('Saving polygons...')
        writer = output.getVectorWriter(fields, QgsWkbTypes.Polygon, vlayer.crs())
        outFeat = QgsFeature()
        total = 50.0 / len(polygons)
        for current, polygon in enumerate(polygons):
            outFeat.setGeometry(QgsGeometry.fromWkt(polygon.wkt))
            if self.getParameterValue(self.GEOMETRY):
                outFeat.setAttributes([None] * fieldsCount + [polygon.area,
                                                              polygon.length])
            writer.addFeature(outFeat)
            progress.setPercentage(50 + int(current * total))
        del writer

    def _cont_to_polys(self, cont_lats, cont_lons, cont_val):
        """
        Take the lat/lon contours, split them into their different segments, and create polygons out of them. Contours
        that stretch from border to border will end up covering large sections of the country. That's okay; we'll take
        care of that later.
        """
        polys = {}

        start_idx = 0
        splits = []
        while True:
            try:
                split_idx = cont_lats[start_idx:].index(99.99)
                splits.append(start_idx + split_idx)
                start_idx += split_idx + 1
            except ValueError:
                break

        splits = [ -1 ] + splits + [ len(cont_lats) + 1 ]
        poly_lats = [ cont_lats[splits[i] + 1:splits[i + 1]] for i in xrange(len(splits) - 1) ]
        poly_lons = [ cont_lons[splits[i] + 1:splits[i + 1]] for i in xrange(len(splits) - 1) ]

        # Intersect with the US boundary shape file.
        for plat, plon in zip(poly_lats, poly_lons):
            cont = LineString(zip(plon, plat))

            if plat[0] != plat[-1] or plon[0] != plon[-1]:
                # If the line is not a closed contour, then it intersects with the edge of the US. Extend
                #   the ends a little bit to make sure it's outside the edge.
                dln = np.diff(plon)
                dlt = np.diff(plat)

                pre = [ (plon[0] - 0.03 * dln[0], plat[0] - 0.03 * dlt[0]) ]
                post = [ (plon[-1] + 0.03 * dln[-1], plat[-1] + 0.03 * dlt[-1]) ]

                cont.coords = pre + list(cont.coords) + post 

            # polygonize() will split the country into two parts: one inside the outlook and one outside.
            #   Construct test_ln that is to the right of (inside) the contour and keep only the polygon
            #   that contains the line
            test_ln = cont.parallel_offset(0.05, 'right')

            polys[cont] = []

            for poly in polygonize(self._conus.boundary.union(cont)):
                if (poly.crosses(test_ln) or poly.contains(test_ln)) and self._conus.contains(poly.buffer(-0.01)):
                    polys[cont].append(poly)

        return polys

  def union(self):
    '''new method using Union function'''
    global polyCount
    inFeat = QgsFeature()

    progress = 0.

    self.emit(SIGNAL('progress'), 0)

    self.parent.layer = getMapLayerByName(self.ui.cmbLayer.currentText())
    provider = self.parent.layer.dataProvider()
    #user can't toggle edit mode of line layer while polygonizing, plugin automatically turn it off 
    QObject.connect(self.parent.layer,SIGNAL("editingStarted()"), self.parent.startEditing)
    allAttrs = provider.attributeIndexes()
    provider.select(allAttrs)

    provider.select()

    step = 45. / self.parent.layer.featureCount()
    allLinesList = []
    allLinesListExtend = allLinesList.extend
    allLinesListAppend = allLinesList.append

    while provider.nextFeature(inFeat):
      geom = inFeat.geometry()
      if geom.isMultipart():
        allLinesListExtend(geom.asMultiPolyline() )
      else:
        allLinesListAppend(geom.asPolyline())

      progress += step
      self.emit(SIGNAL('progress'), progress)

    allLines = MultiLineString(allLinesList)
    allLines = allLines.union(Point(0,0))

    polygons = list(polygonize([allLines]))

    polyCount = len(polygons)
    #if no polygons where created then exit from thread
    if polyCount == 0:
      QObject.disconnect(self.parent.polygonizeThread,SIGNAL("finished()"), self.parent.threadFinished)
      self.emit(SIGNAL('noPolygons'))
      return
    else:
      if self.ui.cbOutput.isChecked():
        self.saveAsFile(polygons, progress)
      else:
        self.saveInMemory(polygons, progress)

    def _geometric_process(self):

        if len(self._entries) != 1:
            merged_union = linemerge(unary_union(self._entries))
        else:
            merged_union = MultiLineString(self._entries)

        self._edges = [ Edge(tuple(g.coords)) for g in merged_union.geoms ]

        edges = map(lambda e: e._geom, self._edges)

        if self._bnode:
            xyset = set([edge.coords[i] for edge in edges for i in (0,-1)])
            self._nodes = map(lambda xy: Node(xy),xyset)

        if self._bface:
            self._faces = [Face(poly.exterior,poly.interiors) for poly in polygonize(edges)]