def AttemptCollapse(name, zone):
  overlaps = False
  for i in range(0,len(zone)):
    for j in range(i+1, len(zone)):
      if zone[i].overlaps(zone[j]):
        overlaps = True
        break
    if overlaps:
      break

  if overlaps:
    print 'Found overlapping zone on %s!' % name
    mp = SMultiPolygon(zone)
    collapsed = cascaded_union(mp)
    if type(collapsed) == SPolygon:
      return [collapsed]
    else:
      print 'Got multipolygon len %d' % len(collapsed.geoms)
      return collapsed.geoms
  else:
    print 'Non-overlapping zone on %s (%d)' % (name, len(zone))
    return zone 

def get_point_in_polygon_score(self, lane_seq: List[int],
                                   xy_seq: np.ndarray, city_name: str,
                                   avm: ArgoverseMap) -> int:
        """Get the number of coordinates that lie insde the lane seq polygon.

        Args:
            lane_seq: Sequence of lane ids
            xy_seq: Trajectory coordinates
            city_name: City name (PITT/MIA)
            avm: Argoverse map_api instance
        Returns:
            point_in_polygon_score: Number of coordinates in the trajectory that lie within the lane sequence

        """
        lane_seq_polygon = cascaded_union([
            Polygon(avm.get_lane_segment_polygon(lane, city_name)).buffer(0)
            for lane in lane_seq
        ])
        point_in_polygon_score = 0
        for xy in xy_seq:
            point_in_polygon_score += lane_seq_polygon.contains(Point(xy))
        return point_in_polygon_score 

def cascaded_union_with_alt(polyalt: AirspaceList) -> AirspaceList:
    altitudes = set(alt for _, *low_up in polyalt for alt in low_up)
    slices = sorted(altitudes)
    if len(slices) == 1 and slices[0] is None:
        simple_union = cascaded_union([p for p, *_ in polyalt])
        return [ExtrudedPolygon(simple_union, float("-inf"), float("inf"))]
    results: List[ExtrudedPolygon] = []
    for low, up in zip(slices, slices[1:]):
        matched_poly = [
            p
            for (p, low_, up_) in polyalt
            if low_ <= low <= up_ and low_ <= up <= up_
        ]
        new_poly = ExtrudedPolygon(cascaded_union(matched_poly), low, up)
        if len(results) > 0 and new_poly.polygon.equals(results[-1].polygon):
            merged = ExtrudedPolygon(new_poly.polygon, results[-1].lower, up)
            results[-1] = merged
        else:
            results.append(new_poly)
    return results


# -- Methods below are placed here because of possible circular imports -- 

def alpha_shape(points, alpha):
    """
    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

    coords = np.array([point.coords[0] for point in points])
    tri = Delaunay(coords)
    triangles = coords[tri.vertices]
    a = ((triangles[:,0,0] - triangles[:,1,0]) ** 2 + (triangles[:,0,1] - triangles[:,1,1]) ** 2) ** 0.5
    b = ((triangles[:,1,0] - triangles[:,2,0]) ** 2 + (triangles[:,1,1] - triangles[:,2,1]) ** 2) ** 0.5
    c = ((triangles[:,2,0] - triangles[:,0,0]) ** 2 + (triangles[:,2,1] - triangles[:,0,1]) ** 2) ** 0.5
    s = ( a + b + c ) / 2.0
    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
    circums = a * b * c / (4.0 * areas)
    filtered = triangles[circums < (1.0 / alpha)]
    edge1 = filtered[:,(0,1)]
    edge2 = filtered[:,(1,2)]
    edge3 = filtered[:,(2,0)]
    edge_points = np.unique(np.concatenate((edge1,edge2,edge3)), axis = 0).tolist()
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return cascaded_union(triangles), edge_points 

def cutpath(self):
        selected = self.get_selected()
        tools = selected[1:]
        toolgeo = cascaded_union([shp.geo for shp in tools])

        target = selected[0]
        if type(target.geo) == Polygon:
            for ring in poly2rings(target.geo):
                self.add_shape(DrawToolShape(ring.difference(toolgeo)))
            self.delete_shape(target)
        elif type(target.geo) == LineString or type(target.geo) == LinearRing:
            self.add_shape(DrawToolShape(target.geo.difference(toolgeo)))
            self.delete_shape(target)
        else:
            self.app.log.warning("Not implemented.")

        self.replot() 

def union(self):
        """
        Makes union of selected polygons. Original polygons
        are deleted.

        :return: None.
        """

        results = cascaded_union([t.geo for t in self.get_selected()])

        # Delete originals.
        for_deletion = [s for s in self.get_selected()]
        for shape in for_deletion:
            self.delete_shape(shape)

        # Selected geometry is now gone!
        self.selected = []

        self.add_shape(DrawToolShape(results))

        self.replot() 

def buffer(self, buf_distance):
        selected = self.get_selected()

        if len(selected) == 0:
            self.app.inform.emit("[warning] Nothing selected for buffering.")
            return

        if not isinstance(buf_distance, float):
            self.app.inform.emit("[warning] Invalid distance for buffering.")
            return

        pre_buffer = cascaded_union([t.geo for t in selected])
        results = pre_buffer.buffer(buf_distance)
        self.add_shape(DrawToolShape(results))

        self.replot() 

def subtract_polygon(self, points):
        """
        Subtract polygon from the given object. This only operates on the paths in the original geometry, i.e. it converts polygons into paths.

        :param points: The vertices of the polygon.
        :return: none
        """
        if self.solid_geometry is None:
            self.solid_geometry = []

        #pathonly should be allways True, otherwise polygons are not subtracted
        flat_geometry = self.flatten(pathonly=True)
        log.debug("%d paths" % len(flat_geometry))
        polygon=Polygon(points)
        toolgeo=cascaded_union(polygon)
        diffs=[]
        for target in flat_geometry:
            if type(target) == LineString or type(target) == LinearRing:
                diffs.append(target.difference(toolgeo))
            else:
                log.warning("Not implemented.")
        self.solid_geometry=cascaded_union(diffs) 

def export_svg(self, scale_factor=0.00):
        """
        Exports the Gemoetry Object as a SVG Element

        :return: SVG Element
        """
        # Make sure we see a Shapely Geometry class and not a list
        geom = cascaded_union(self.flatten())

        # scale_factor is a multiplication factor for the SVG stroke-width used within shapely's svg export

        # If 0 or less which is invalid then default to 0.05
        # This value appears to work for zooming, and getting the output svg line width
        # to match that viewed on screen with FlatCam
        if scale_factor <= 0:
            scale_factor = 0.05

        # Convert to a SVG
        svg_elem = geom.svg(scale_factor=scale_factor)
        return svg_elem 

def _get_merged_polygon(self, cluster):
        """Merge polygons using shapely (Internal).

        Given a single cluster from _get_merge_clusters_from_df(), This creates
        and merges polygons into a single cascaded union. It first dilates the
        polygons by buffer_size pixels to make them overlap, merges them,
        then erodes back by buffer_size to get the merged polygon.

        """
        buffer_size = self.merge_thresh + 3
        nest_polygons = []
        for nestinfo in cluster:
            nest = dict(self.edge_contours[nestinfo['roiname']].loc[
                nestinfo['nid'], :])
            roitop = self.roiinfos[nestinfo['roiname']]['top']
            roileft = self.roiinfos[nestinfo['roiname']]['left']
            coords = _parse_annot_coords(
                nest, x_offset=roileft, y_offset=roitop)
            nest_polygons.append(Polygon(coords).buffer(buffer_size))
        merged_polygon = cascaded_union(nest_polygons).buffer(-buffer_size)
        return merged_polygon

    # %% ===================================================================== 

def LinesToPolyHoles(lines):
    # get all of the polys, both outer and holes
    # if there are holes, you'll get them double:
    # once as themselves and again as holes in a boundary
    polys = list(polygonize(lines))

    # merge to get just the outer
    bounds =  cascaded_union(polys)

    if not iterable(bounds):
        bounds = [bounds]

    retval = []
    for bound in bounds:
        for poly in polys:
            if Polygon(poly.exterior).almost_equals(bound):
                retval.append(poly)

    return retval 

def offset(self, dist, side):
        if dist == 0:
            return _SidedPolygon(self.poly, self.level)
        if dist < 0:
            side = not side
            dist = abs(dist)
        if (side == c.OUTSIDE and self.isFeature) or (side == c.INSIDE and not self.isFeature):
            return _SidedPolygon(self.poly.buffer(dist), self.level)
        try:
            buffPoly = self.poly.exterior.buffer(dist)
            if len(buffPoly.interiors) > 1:
                inPoly = cascaded_union([Polygon(i) for i in buffPoly.interiors])            
            else:
                inPoly = Polygon(buffPoly.interiors[0])
            return _SidedPolygon(inPoly, self.level)
        except Exception:
            return None 

def _area_at_zoom(self, zoom):
        if zoom not in self._cache_area_at_zoom:
            # use union of all input items and, if available, intersect with
            # init_bounds
            if "input" in self._params_at_zoom[zoom]:
                input_union = cascaded_union([
                    self.input[get_hash(v)].bbox(self.process_pyramid.crs)
                    for k, v in _flatten_tree(self._params_at_zoom[zoom]["input"])
                    if v is not None
                ])
                self._cache_area_at_zoom[zoom] = input_union.intersection(
                    box(*self.init_bounds)
                ) if self.init_bounds else input_union
            # if no input items are available, just use init_bounds
            else:
                self._cache_area_at_zoom[zoom] = box(*self.init_bounds)
        return self._cache_area_at_zoom[zoom] 

def get_names_rings(area):
    names = []
    all_rings =[]

    for idx, a in enumerate(area):
        rings = [LinearRing(p) for p in a['polygons']]
        if len(rings) > 1:
            u = cascaded_union(rings)
        else:
            u = rings[0]
        all_rings.append(u.envelope)
        if a.get('name'):
            names.append(a.get('name'))

    names = sorted(names)
    
    return names, all_rings 

def ReadShapeFile(ShapeFile):

    """
    Open shapefile and create Polygon dictionary
    returns dictionary of shapely Polygons

    MDH

    """

    #open shapefile and read shapes
    Shapes = fiona.open(ShapeFile)

    # get the input coordinate system
    Input_CRS = Proj(Shapes.crs)

    # Create a dictionary of shapely polygons
    PolygonDict = {}

    # loop through shapes and add to dictionary
    for Feature in Shapes:
        if Feature['geometry']['type'] == 'MultiPolygon':
            Shape = MultiPolygon(shape(Feature['geometry']))
            Value = float(Feature['properties']['ID'])
        elif Feature['geometry']['type'] == 'Polygon':
            Shape = Polygon(shape(Feature['geometry']))
            Value = float(Feature['properties']['ID'])

        #check for multipolygons
        if Value in PolygonDict:
            Polygons = [Shape, PolygonDict[Value]]
            Shape = cascaded_union(Polygons)
        #update dictionary
        PolygonDict[Value] = Shape

    return PolygonDict, Input_CRS 

def _polytree_to_shapely(tree):
    polygons, children = _polytree_node_to_shapely(tree)

    # expect no left over children - should all be incorporated into polygons
    # by the time recursion returns to the root.
    assert len(children) == 0

    union = cascaded_union(polygons)
    assert union.is_valid
    return union 

def isPpaWithinServiceArea(pal_records, ppa_zone_geometry):
  """Check if the ppa zone geometry with in service area then return True.

  Checks the ppa zone geometry's boundary and interior intersect only with the
  interior of the service area (not its boundary or exterior).

  Args:
    pal_records: A list of pal records to compute service area based on
      census_tracts.
    ppa_zone_geometry: A PPA polygon dictionary in GeoJSON format.

  Returns:
    A value is the boolean with the value as True if the ppa zone geometry's
      boundary and interior intersect with in the interior of the service
      area otherwise value as false.

  """
  # Get the census tract for each pal record and convert it to Shapely
  # geometry.
  census_tracts_for_pal = [
      utils.ToShapely(drive.census_tract_driver.GetCensusTract(
          pal['license']['licenseAreaIdentifier'])
                      ['features'][0]['geometry']) for pal in pal_records]
  pal_service_area = ops.cascaded_union(census_tracts_for_pal)

  # Convert GeoJSON dictionary to Shapely object.
  ppa_zone_shapely_geometry = utils.ToShapely(ppa_zone_geometry)

  return ppa_zone_shapely_geometry.buffer(-1e-6).within(pal_service_area) 

def _ClipPpaByCensusTract(contour_union, pal_records):
  """ Clip a PPA 'contour_union' zone (shapely.MultiPolygon)
  with the census tracts defined by a sequence of 'pal_records'."""

  # Get the Census Tract for Each Pal Record and Convert it to Shapely Geometry.
  census_tracts_for_pal = [sgeo.shape(
      drive.census_tract_driver.GetCensusTract(pal['license']['licenseAreaIdentifier'])
      ['features'][0]['geometry']).buffer(0) for pal in pal_records]
  census_tracts_union = ops.cascaded_union(census_tracts_for_pal)
  return contour_union.intersection(census_tracts_union) 

def PpaCreationModel(devices, pal_records):
  """Creates a PPA Polygon based on the PAL Records and Device Information
  Args:
    devices: A list of CBSD records (schema |CbsdRecordData|).
    pal_records: A list of pal records.
  Returns:
    The PPA polygon in GeoJSON format (string).
  """
  # Validation for Inputs
  for device in devices:
    logging.info('Validating device', device)
    util.assertContainsRequiredFields("RegistrationRequest.schema.json", device)
  for pal_rec in pal_records:
    logging.info('Validating pal_rec', pal_rec)
    util.assertContainsRequiredFields("PalRecord.schema.json", pal_rec)

  # Create Contour for each CBSD
  pool = mpool.Pool()
  device_polygon = pool.map(_GetPolygon, devices)

  # Create Union of all the CBSD Contours and Check for hole
  # after Census Tract Clipping
  contour_union = ops.cascaded_union(device_polygon)
  logging.info('contour_union = %s', contour_union)

  ppa_without_small_holes = utils.PolyWithoutSmallHoles(contour_union)
  logging.info('ppa_without_small_holes = %s', ppa_without_small_holes)

  ppa_polygon = _ClipPpaByCensusTract(ppa_without_small_holes, pal_records)
  logging.info('contour_union clipped by census tracts: %s', ppa_polygon)
  if ppa_polygon.is_empty:
    raise Exception("Empty Polygon is generated, please check the inputs.")
  if ppa_polygon.geom_type == "MultiPolygon":
    raise Exception("Multi Polygon is not supported, please check the inputs.")

  # Convert Shapely Object to GeoJSON geometry string
  return utils.ToGeoJson(ppa_polygon) 

def computeUnionArea(boxes):
    boxes = [geo.box(bb[0],bb[1],bb[0]+bb[2], bb[1]+bb[3]) for bb in boxes]
    return cascaded_union(boxes).area 

def fill_waveguide_with_holes_in_honeycomb_lattice(waveguide, spacing, padding, hole_radius):
    """
    Fills a given waveguide with holes which are arranged in a honeycomb structure
    This can be used for generating vortex traps like presented in https://doi.org/10.1103/PhysRevApplied.11.064053

    :param waveguide: Waveguide to be filled
    :param spacing: Spacing between the holes
    :param padding: Minimum distance from the edge of the waveguide to the holes
    :param hole_radius: Radius of the holes
    :return: Shapely object, which describes the holes
    """
    center_coordinates = LineString(waveguide.center_coordinates)
    outline = prep(waveguide.get_shapely_outline().buffer(hole_radius).buffer(-padding - 2 * hole_radius))
    area_for_holes = prep(waveguide.get_shapely_object().buffer(hole_radius).buffer(-padding - 2 * hole_radius))
    circles = []

    offset = 0
    new_circles = True
    while new_circles:
        new_circles = False
        for i, pos in enumerate(np.arange(padding, center_coordinates.length - padding, np.sqrt(3 / 4) * spacing)):
            xy = np.array(center_coordinates.interpolate(pos))
            diff = np.array(center_coordinates.interpolate(pos + padding / 2)) - np.array(
                center_coordinates.interpolate(pos - padding / 2))
            d1 = np.array((-diff[1], diff[0])) / np.linalg.norm(diff)
            for direction in [-1, 1]:
                point = Point(xy + direction * (offset + spacing / 2 * (i % 2)) * d1)
                if outline.contains(point):
                    new_circles = True
                if area_for_holes.contains(point):
                    circles.append(point.buffer(hole_radius))
        offset += spacing
    return cascaded_union(circles) 

def get_shapely_object(self):
        # Start with the "rays"
        # noinspection PyTypeChecker
        ray_angles = np.linspace(0, np.pi / 2, 8) + np.pi / 2
        outer_ray_positions = (np.array([np.cos(ray_angles), np.sin(ray_angles)]) * self.height).T + (self.height, 0)
        inner_ray_positions = (np.array([np.cos(ray_angles), np.sin(ray_angles)])
                               * self.height * self.min_radius_fraction).T + (self.height, 0)

        polygons = list()
        for outer, inner in zip(outer_ray_positions.reshape(-1, 2, 2), inner_ray_positions.reshape(-1, 2, 2)):
            polygons.append(Polygon([outer[0], outer[1], inner[1], inner[0]]))
            pass

        # Draw the letters
        d = self.height * 0.077

        k_upper_branch = rotate(box(0, -d, sqrt(2) * self.height / 2 + sqrt(2) * d, d), 45, origin=(0, 0))
        k_lower_branch = scale(k_upper_branch, yfact=-1., origin=(0, 0))
        k_uncut = k_upper_branch.union(k_lower_branch)
        k_unscaled = k_uncut.intersection(box(0, -self.height / 2., self.height / 2. + sqrt(2) * d, self.height / 2.))
        k = scale(k_unscaled, 0.8, origin=(0, 0))
        polygons.append(translate(k, self.height * 1.05, self.height / 2.))

        i = box(0, 0, 2 * d, self.height)
        polygons.append(translate(i, self.height * 1.6))

        t_overlap = 2
        t = box(-d, 0, d, self.height).union(
            box(-d * (1 + t_overlap), self.height - 2 * d, d * (1 + t_overlap), self.height))
        polygons.append(translate(t, self.height * 2.05))

        logo = cascaded_union(polygons)

        return translate(logo, *self.origin) 

def _polytree_to_shapely(tree):
    polygons, children = _polytree_node_to_shapely(tree)

    # expect no left over children - should all be incorporated into polygons
    # by the time recursion returns to the root.
    assert len(children) == 0

    union = cascaded_union(polygons)
    assert union.is_valid
    return union 

def flatten(self) -> Polygon:
        """Returns the 2D footprint of the airspace."""
        return cascaded_union([p.polygon for p in self]) 

def shape(self):
        # filter out the contour, helps for useful display
        # list(self.osm_request()),
        return self.osm_request().shape
        return cascaded_union(
            list(
                shape
                for dic, shape in self.osm_request().ways.values()
                if dic["tags"]["aeroway"] != "aerodrome"
            )
        ) 

def subtract(self):
        selected = self.get_selected()
        tools = selected[1:]
        toolgeo = cascaded_union([shp.geo for shp in tools])
        result = selected[0].geo.difference(toolgeo)

        self.delete_shape(selected[0])
        self.add_shape(DrawToolShape(result))

        self.replot() 

def import_svg(self, filename, flip=True):
        """
        Imports shapes from an SVG file into the object's geometry.

        :param filename: Path to the SVG file.
        :type filename: str
        :return: None
        """

        # Parse into list of shapely objects
        svg_tree = ET.parse(filename)
        svg_root = svg_tree.getroot()

        # Change origin to bottom left
        # h = float(svg_root.get('height'))
        # w = float(svg_root.get('width'))
        h = svgparselength(svg_root.get('height'))[0]  # TODO: No units support yet
        geos = getsvggeo(svg_root)

        if flip:
            geos = [translate(scale(g, 1.0, -1.0, origin=(0, 0)), yoff=h) for g in geos]

        # Add to object
        if self.solid_geometry is None:
            self.solid_geometry = []

        if type(self.solid_geometry) is list:
            self.solid_geometry.append(cascaded_union(geos))
        else:  # It's shapely geometry
            self.solid_geometry = cascaded_union([self.solid_geometry,
                                                  cascaded_union(geos)])

        return 

def union(self):
        """
        Runs a cascaded union on the list of objects in
        solid_geometry.

        :return: None
        """
        self.solid_geometry = [cascaded_union(self.solid_geometry)] 

def make_moire(mods):
        """
        Note: Specs indicate that rotation is only allowed if the center
        (x, y) == (0, 0). I will tolerate breaking this rule.

        :param mods: (x-center, y-center, outer_dia_outer_ring, ring thickness,
            gap, max_rings, crosshair_thickness, crosshair_len, rotation
            angle around origin in degrees)
        :return:
        """

        x, y, dia, thickness, gap, nrings, cross_th, cross_len, angle = ApertureMacro.default2zero(9, mods)

        r = dia/2 - thickness/2
        result = Point((x, y)).buffer(r).exterior.buffer(thickness/2.0)
        ring = Point((x, y)).buffer(r).exterior.buffer(thickness/2.0)  # Need a copy!

        i = 1  # Number of rings created so far

        ## If the ring does not have an interior it means that it is
        ## a disk. Then stop.
        while len(ring.interiors) > 0 and i < nrings:
            r -= thickness + gap
            if r <= 0:
                break
            ring = Point((x, y)).buffer(r).exterior.buffer(thickness/2.0)
            result = cascaded_union([result, ring])
            i += 1

        ## Crosshair
        hor = LineString([(x - cross_len, y), (x + cross_len, y)]).buffer(cross_th/2.0, cap_style=2)
        ver = LineString([(x, y-cross_len), (x, y + cross_len)]).buffer(cross_th/2.0, cap_style=2)
        result = cascaded_union([result, hor, ver])

        return {"pol": 1, "geometry": result} 

def create_geometry(self):
        # TODO: This takes forever. Too much data?
        self.solid_geometry = cascaded_union([geo['geom'] for geo in self.gcode_parsed])