Python shapely.geometry 模块，Polygon() 实例源码

def test_intersection_com_mock_2(self):
        ls = LineString([(1, 1, 9.48024060e+08), (2, 2, 9.49363260e+08),
                         (3, 1, 9.51868860e+08)])

        poly = Polygon([(1, 1), (1, 3), (4, 3), (4, 1), (1, 1)])

        self.traj2.intersection_shapely = MagicMock(return_value=ls)
        response = self.traj2.intersection_shapely(poly)

        ls = np.array(ls)

        trajMock = self.traj2.to_Trajectory(response)
        traj = Trajectory(ls[:, 0], ls[:, 1], ls[:, 2])

        assert (np.array_equal(trajMock.getX(), traj.getX()))
        assert (np.array_equal(trajMock.getY(), traj.getY()))
        assert (np.array_equal(trajMock.getTime(), traj.getTime()))

def generate_voronoi(geoseries_polygons):
    """Generate Voronoi polygons from polygon edges
    :param geoseries_polygons: GeoSeries of raw polygons
    :return:
    """
    edges = geoseries_polygons.unary_union.boundary
    pnts = points_along_boundaries(edges, 0.75)
    cent = pnts.centroid
    tpnts = translate(pnts, -cent.x, -cent.y)
    vor = Voronoi(tpnts)
    polys = []
    for region in vor.regions:
        if len(region) > 0 and all([i > 0 for i in region]):
            polys.append(Polygon([vor.vertices[i] for i in region]))
    gs_vor = geopandas.GeoSeries(polys)
    t_gs_vor = gs_vor.translate(cent.x, cent.y)
    t_gs_vor.crs = geoseries_polygons.crs
    return t_gs_vor, pnts

def bounded_segments(lines, bounding_box, cut_segment=True):
    """
    Extract the bounded segments from a list of lines
    :param lines: a list of LineString
    :param bounding_box: the bounding coordinates in (minx, miny, maxx, maxy)
           or Polygon instance
    :return: a list of bounded segments
    """
    if isinstance(bounding_box, Polygon):
        bbox = bounding_box
    else:
        bbox = box(bounding_box[0], bounding_box[1],
                   bounding_box[2], bounding_box[3])
    segments = []
    for line in lines:
        if line.intersects(bbox):
            if cut_segment:
                segments.append(line.intersection(bbox))
            else:
                segments.append(line)
    return segments

def subgraph_within_box(self, bounding_box):
        """
        Extract a subgraph bounded by a box.
        :param bounding_box: the bounding coordinates in
            (minx, miny, maxx, maxy) or a Polygon instance
        :return: a subgraph of nx.Graph
        """
        if isinstance(bounding_box, Polygon):
            bbox = bounding_box
        else:
            bbox = box(bounding_box[0], bounding_box[1],
                       bounding_box[2], bounding_box[3])
        nbunch = set()
        for edge in self.graph.edges():
            s, e = edge
            if bbox.intersects(LineString([self.node_xy[s], self.node_xy[e]])):
                nbunch.add(s)
                nbunch.add(e)
        return self.graph.subgraph(nbunch)

def polygonal_max(self, geometry, data_type):
        """Finds the max value for each band that is contained within the given geometry.

        Args:
            geometry (shapely.geometry.Polygon or shapely.geometry.MultiPolygon or bytes): A
                Shapely ``Polygon`` or ``MultiPolygon`` that represents the area where the summary
                should be computed; or a WKB representation of the geometry.
            data_type (type): The type of the values within the rasters. Can either be int or
                float.

        Returns:
            [int] or [float] depending on ``data_type``.

        Raises:
            TypeError: If ``data_type`` is not an int or float.
        """

        if data_type is int:
            return self._process_polygonal_summary(geometry, self.srdd.polygonalMax)
        elif data_type is float:
            return self._process_polygonal_summary(geometry, self.srdd.polygonalMaxDouble)
        else:
            raise TypeError("data_type must be either int or float.")

def polygonal_sum(self, geometry, data_type):
        """Finds the sum of all of the values in each band that are contained within the given geometry.

        Args:
            geometry (shapely.geometry.Polygon or shapely.geometry.MultiPolygon or bytes): A
                Shapely ``Polygon`` or ``MultiPolygon`` that represents the area where the summary
                should be computed; or a WKB representation of the geometry.
            data_type (type): The type of the values within the rasters. Can either be int or
                float.

        Returns:
            [int] or [float] depending on ``data_type``.

        Raises:
            TypeError: If ``data_type`` is not an int or float.
        """

        if data_type is int:
            return self._process_polygonal_summary(geometry, self.srdd.polygonalSum)
        elif data_type is float:
            return self._process_polygonal_summary(geometry, self.srdd.polygonalSumDouble)
        else:
            raise TypeError("data_type must be either int or float.")

def polygonize_cv(mask, epsilon=1., min_area=10.):
    contours, hierarchy = cv2.findContours(mask, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_TC89_KCOS)
    # create approximate contours to have reasonable submission size
    approx_contours = [cv2.approxPolyDP(cnt, epsilon, True)
                       for cnt in contours]
    approx_contours = contours
    if not contours:
        return MultiPolygon()
    # now messy stuff to associate parent and child contours
    cnt_children = defaultdict(list)
    child_contours = set()
    assert hierarchy.shape[0] == 1
    # http://docs.opencv.org/3.1.0/d9/d8b/tutorial_py_contours_hierarchy.html
    for idx, (_, _, _, parent_idx) in enumerate(hierarchy[0]):
        if parent_idx != -1:
            child_contours.add(idx)
            cnt_children[parent_idx].append(approx_contours[idx])
    # create actual polygons filtering by area (removes artifacts)
    all_polygons = []
    for idx, cnt in enumerate(approx_contours):
        if idx not in child_contours and cv2.contourArea(cnt) >= min_area:
            assert cnt.shape[1] == 1
            poly = Polygon(
                shell=cnt[:, 0, :],
                holes=[c[:, 0, :] for c in cnt_children.get(idx, [])
                       if cv2.contourArea(c) >= min_area])
            all_polygons.append(poly)
    # approximating polygons might have created invalid ones, fix them
    all_polygons = MultiPolygon(all_polygons)
    if not all_polygons.is_valid:
        all_polygons = all_polygons.buffer(0)
        # Sometimes buffer() converts a simple Multipolygon to just a Polygon,
        # need to keep it a Multi throughout
        if all_polygons.type == 'Polygon':
            all_polygons = MultiPolygon([all_polygons])
    return all_polygons

def test_intersection_com_mock(self):
        ls = LineString([(1.5, 1, 9.48024060e+08), (2, 2, 9.49363260e+08),
                         (3, 2, 9.51868860e+08), (4, 3, 9.53208060e+08)])

        poly = Polygon([(1, 1), (1, 3), (4, 3), (4, 1), (1, 1)])

        self.traj.intersection_shapely = MagicMock(return_value=ls)
        response = self.traj.intersection_shapely(poly)

        ls = np.array(ls)
        trajMock = self.traj.to_Trajectory(response)
        traj = Trajectory(ls[:, 0], ls[:, 1], ls[:, 2])

        assert (np.array_equal(trajMock.getX(), traj.getX()))
        assert (np.array_equal(trajMock.getY(), traj.getY()))
        assert (np.array_equal(trajMock.getTime(), traj.getTime()))

def _intersects(self, context, coord, event):
        t = time.time()
        c0 = self._position_3d_from_coord(context, coord)
        c1 = self._position_3d_from_coord(context, (coord[0], event.mouse_region_y))
        c2 = self._position_3d_from_coord(context, (event.mouse_region_x, event.mouse_region_y))
        c3 = self._position_3d_from_coord(context, (event.mouse_region_x, coord[1]))
        poly = ShapelyPolygon([c0, c1, c2, c3])
        prepared_poly = shapely.prepared.prep(poly)
        count, gids = self.tree.intersects(poly)
        if event.ctrl:
            selection = [i for i in gids if prepared_poly.contains(self.geoms[i])]
        else:
            selection = [i for i in gids if prepared_poly.intersects(self.geoms[i])]
        print("Selectable._intersects() :%.2f seconds" % (time.time() - t))
        if event.shift:
            self._unselect(selection)
        else:
            self._select(selection)
        self._draw(context)

def _process_element(self, element):
        if element.interface.datatype == 'geodataframe':
            geoms = element.split(datatype='geom')
            projected = [self.p.projection.project_geometry(geom, element.crs)
                         for geom in geoms]
            new_data = element.data.copy()
            new_data['geometry'] = projected
            return element.clone(new_data, crs=self.p.projection)

        geom_type = Polygon if isinstance(element, Polygons) else LineString
        xdim, ydim = element.kdims[:2]
        projected = []
        for geom in element.split(datatype='columns'):
            xs, ys = geom[xdim.name], geom[ydim.name]
            path = geom_type(np.column_stack([xs, ys]))
            proj = self.p.projection.project_geometry(path, element.crs)
            proj_arr = geom_to_array(proj)
            geom[xdim.name] = proj_arr[:, 0]
            geom[ydim.name] = proj_arr[:, 1]
            projected.append(geom)
        return element.clone(projected, crs=self.p.projection)

def get_contour_mask(dd, id, dosegridpoints, contour):
    """Get the mask for the contour with respect to the dose plane."""
    doselut = dd['lut']

    c = matplotlib.path.Path(list(contour))

    # def inpolygon(polygon, xp, yp):
    #     return np.array(
    #         [Point(x, y).intersects(polygon) for x, y in zip(xp, yp)],
    #         dtype=np.bool)

    # p = Polygon(contour)
    # x, y = np.meshgrid(np.array(dd['lut'][0]), np.array(dd['lut'][1]))
    # mask = inpolygon(p, x.ravel(), y.ravel())
    # return mask.reshape((len(doselut[1]), len(doselut[0])))

    grid = c.contains_points(dosegridpoints)
    grid = grid.reshape((len(doselut[1]), len(doselut[0])))

    return grid

def Import_reticule(IDT_data,reticule_size,reticule_filename):
    all_points = SVGT.read_path_from_svg(reticule_filename)
    reticule =[]
    for points in all_points:
        reticule.append(shapely_geom.Polygon(points.T))
    reticule = shapely_geom.MultiPolygon(reticule)
    outbox = reticule.bounds
    dx = outbox[2]-outbox[0]
    dy = outbox[3]-outbox[1]
    d = max([dx,dy])
    x0 = outbox[0]+dx/2
    y0 = outbox[1]+dy/2
    factor = 2*reticule_size/d
    reticule = shapely_affinity.translate(reticule, xoff=-x0, yoff=-y0) 
    reticule = shapely_affinity.scale(reticule, xfact = factor, yfact= factor, origin=(0,0,0))
    IDT_data['reticule'] = reticule

def Import_SVG_route(IDT_data):
    route_filename = IDT_data['svg_route']
    all_points = SVGT.read_path_from_svg(route_filename)
    route =[]
    for points in all_points:
        route.append(shapely_geom.Polygon(points.T))
    route = shapely_geom.MultiPolygon(route)
    #outbox = route.bounds
    #dx = outbox[2]-outbox[0]
    #dy = outbox[3]-outbox[1]
    #x0 = outbox[0]+dx/2
    #y0 = outbox[1]+dy/2
    x0svg = 4000
    y0svg = 4000
    factor = 1e-5;
    route = shapely_affinity.translate(route, xoff=-x0svg, yoff=-y0svg) 
    route = shapely_affinity.scale(route, xfact = factor, yfact= factor, origin=(0,0,0))
    IDT_data['route'] = route

def test_slice_linearring(self):
        shell = geometry.polygon.LinearRing([(0.0, 0.0), (70.0, 120.0),
                                             (140.0, 0.0), (0.0, 0.0)])
        holes = [geometry.polygon.LinearRing([(60.0, 80.0), (80.0, 80.0),
                                              (70.0, 60.0), (60.0, 80.0)]),
                 geometry.polygon.LinearRing([(30.0, 10.0), (50.0, 10.0),
                                              (40.0, 30.0), (30.0, 10.0)]),
                 geometry.polygon.LinearRing([(90.0, 10), (110.0, 10.0),
                                              (100.0, 30.0), (90.0, 10.0)])]
        g = geometry.Polygon(shell, holes)
        self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[1:], holes[1:])]))
        self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[:-1], holes[:-1])]))
        self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[::-1], holes[::-1])]))
        self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[::2], holes[::2])]))
        self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[:3], holes[:3])]))
        self.assertTrue(g.interiors[3:] == holes[3:] == [])

def get_ways_data(elements, coords):
    ways_data = {}
    for element in elements:
        # Only process ways
        if element.get('type') != 'way':
            continue
        # Only process ways with 3 or more nodes, otherwise
        # Shapely will complain.
        nodes = element.get('nodes')
        if len(nodes) < 3:
            continue
        exterior = [(coords[node].get('lat'), coords[node].get('lon')) \
            for node in nodes]
        # Build the polygon and compute its bbox and centroid
        way_polygon = Polygon(exterior)
        ways_data[element.get('id')] = {
            'bounds': way_polygon.bounds,
            'lat': way_polygon.centroid.x,
            'lon': way_polygon.centroid.y}

    # Done
    return ways_data

def createCirclePolygon(h, k, r, dx):
    """Create shapely polygon of a circle.

    usage: p = createCirclePolygon(h, k, r, dx)

    Args:
        h: x coordinate of center.
        k: y coordinate of center.
        r: radius of circle.
        dx: approximate distance between points * 10.

    Returns:
        Tuple (x, y) of numpy arrays of x and y coordinates of points.
    """
    D = 10.0
    theta = 2 * np.arccos((r - (dx / D)) / r)
    npoints = int(360.0 / theta)
    x, y = getPointsInCircum(r, n=npoints, h=h, k=k)
    p = Polygon(list(zip(x, y)))
    return p

def getClassBalance(pshapes, bounds, proj):
    """
    Get native class balance of projected shapes, assuming a rectangular
    bounding box.

    Args:
        pshapes: Sequence of projected shapely shapes.
        bounds: Desired bounding box, in decimal degrees.
        proj: PyProj object defining orthographic projection of shapes.

    Returns:
        Float fraction of hazard polygons (area of hazard polygons/total
        area of bbox).
    """

    xmin, ymin, xmax, ymax = bounds
    bpoly = Polygon([(xmin, ymax),
                     (xmax, ymax),
                     (xmax, ymin),
                     (xmin, ymin)])
    project = partial(
        pyproj.transform,
        pyproj.Proj(proj='latlong', datum='WGS84'),
        proj)
    bpolyproj = transform(project, bpoly)
    totalarea = bpolyproj.area
    polyarea = 0
    for pshape in pshapes:
        polyarea += pshape.area

    return polyarea / totalarea

def getProjectedShapes(shapes, xmin, xmax, ymin, ymax):
    """
    Take a sequence of geographic shapes and project them to a bounds-centered
    orthographic projection.

    Args:
        shapes: Sequence of shapes, as read in by fiona.collection().
        xmin: Eastern boundary of all shapes.
        xmax: Western boundary of all shapes.
        ymin: Southern boundary of all shapes.
        ymax: Northern boundary of all shapes.

    Returns:
       Tuple of
           - Input sequence of shapes, projected to orthographic
           - PyProj projection object used to transform input shapes
    """
    latmiddle = ymin + (ymax - ymin) / 2.0
    lonmiddle = xmin + (xmax - xmin) / 2.0
    projstr = ('+proj=ortho +datum=WGS84 +lat_0=%.4f +lon_0=%.4f '
               '+x_0=0.0 +y_0=0.0' % (latmiddle, lonmiddle))
    proj = pyproj.Proj(projparams=projstr)
    project = partial(
        pyproj.transform,
        pyproj.Proj(proj='latlong', datum='WGS84'),
        proj)

    pshapes = []
    for tshape in shapes:
        if tshape['geometry']['type'] == 'Polygon':
            pshapegeo = shape(tshape['geometry'])
        else:
            pshapegeo = shape(tshape['geometry'])
        pshape = transform(project, pshapegeo)
        pshapes.append(pshape)  # assuming here that these are simple polygons

    return (pshapes, proj)

def concurrentRequest(self):
        # ???
        envelope=polygon_target.envelope
        bounds=list(envelope.bounds)

        # ????
        bounds[0] -= 0.02
        parts = 4
        # ?????4?????16???
        boundsList = GeoUtil().getBoundsList(bounds, parts)

        threads = []
        for index in range(0, len(boundsList), 1):
            print 'current bounds ...%s ' % index
            subBounds = boundsList[index]

            # ?extent???polygon
            coords=GeoUtil().getPolygonByExtent(subBounds)
            coords=tuple(coords)
            isIntersects=Polygon((coords)).intersects(polygon_target)
            if isIntersects:
                threads.append(gevent.spawn(self.fetchPlace, index, subBounds))
        gevent.joinall(threads)

def mutiSearchPlace(self):
        envelope=polygon_target.envelope
        bounds=list(envelope.bounds)
        # ????
        bounds[0] -= 0.02
        parts = 50
        # ?????4?????16???
        boundsList = GeoUtil().getBoundsList(bounds, parts)
        threads = []
        # ???????????16?????????????
        for index in range(0, len(boundsList)/16+1, 1):
            for threadIndex in range(index*16,(index+1)*16):
                if threadIndex < len(boundsList):
                    print 'current bounds ...%s ' % threadIndex
                    subBounds = boundsList[threadIndex]
                    # ?extent???polygon
                    coords=GeoUtil().getPolygonByExtent(subBounds)
                    coords=tuple(coords)
                    isIntersects=Polygon((coords)).intersects(polygon_target)
                    if isIntersects:
                        threads.append(gevent.spawn(self.fetchPlaceDetail, threadIndex%16, subBounds))
            gevent.joinall(threads)

def error_center_line(self, image, npoints = all, order = 1, overlap = True):
    """Error between worm center line and image

    Arguments:
      image (array): gray scale image of worm
    """
    import shapely.geometry as geo

    if overlap:
      xyl, xyr, cl = self.sides(npoints = npoints);
      w = np.ones(npoints);
      #plt.figure(141); plt.clf();
      worm = geo.Polygon();
      for i in range(npoints-1):
        poly = geo.Polygon(np.array([xyl[i,:], xyr[i,:], xyr[i+1,:], xyl[i+1,:]]));
        ovl = worm.intersection(poly).area;
        tot = poly.area;
        w[i+1] = 1 - ovl / tot;
        worm = worm.union(poly);

      #print w
      return np.sum(w * nd.map_coordinates(image.T, cl.T, order = order))
    else:
      cl = self.center_line(npoints = npoints);
      return np.sum(nd.map_coordinates(image.T, cl.T, order = order))

def contours_from_shape_discrete(left, right):
  """Convert the worm shape to contours

  Arguments:
    left, right (nx2 arrays): left and right side of the worm

  Returns
    nx2: contours of the worm
  """
  poly = geom.Polygon(np.vstack([left, right[::-1,:]]));
  poly = poly.buffer(0)
  bdr = poly.boundary;
  if isinstance(bdr, geom.multilinestring.MultiLineString):
    cts = [];
    for b in bdr:
      x,y = b.xy;
      cts.append(np.vstack([x,y]).T);
  else: # no self intersections
    x,y = bdr.xy;
    cts = np.vstack([x,y]).T;

  return tuple(cts)

def concave_hull(points, alpha, delunay_args=None):
    """Computes the concave hull (alpha-shape) of a set of points.
    """
    delunay_args = delunay_args or {
        'furthest_site': False,
        'incremental': False,
        'qhull_options': None
    }
    triangulation = Delaunay(np.array(points))
    alpha_complex = get_alpha_complex(alpha, points, triangulation.simplices)
    X, Y = [], []
    for s in triangulation.simplices:
        X.append([points[s[k]][0] for k in [0, 1, 2, 0]])
        Y.append([points[s[k]][1] for k in [0, 1, 2, 0]])
    poly = Polygon(list(zip(X[0], Y[0])))
    for i in range(1, len(X)):
        poly = poly.union(Polygon(list(zip(X[i], Y[i]))))
    return poly

def _intersects(self, context, coord, event):
        t = time.time()
        c0 = self._position_3d_from_coord(context, coord)
        c1 = self._position_3d_from_coord(context, (coord[0], event.mouse_region_y))
        c2 = self._position_3d_from_coord(context, (event.mouse_region_x, event.mouse_region_y))
        c3 = self._position_3d_from_coord(context, (event.mouse_region_x, coord[1]))
        poly = ShapelyPolygon([c0, c1, c2, c3])
        prepared_poly = shapely.prepared.prep(poly)
        count, gids = self.tree.intersects(poly)
        if event.ctrl:
            selection = [i for i in gids if prepared_poly.contains(self.geoms[i])]
        else:
            selection = [i for i in gids if prepared_poly.intersects(self.geoms[i])]
        print("Selectable._intersects() :%.2f seconds" % (time.time() - t))
        if event.shift:
            self._unselect(selection)
        else:
            self._select(selection)
        self._draw(context)

def buildings_from_polygon(polygon, retain_invalid=False):
    """
    Get building footprints within some polygon.

    Parameters
    ----------
    polygon : Polygon
    retain_invalid : bool
        if False discard any building footprints with an invalid geometry

    Returns
    -------
    GeoDataFrame
    """

    return create_buildings_gdf(polygon=polygon, retain_invalid=retain_invalid)

def removeIgnoredPoints(gtFramesAll,prFramesAll):

    imgidxs = []
    for imgidx in range(len(gtFramesAll)):
        if ("ignore_regions" in gtFramesAll[imgidx].keys() and
            len(gtFramesAll[imgidx]["ignore_regions"]) > 0):
            regions = gtFramesAll[imgidx]["ignore_regions"]
            polyList = []
            for ridx in range(len(regions)):
                points = regions[ridx]["point"]
                pointList = []
                for pidx in range(len(points)):
                    pt = geometry.Point(points[pidx]["x"][0], points[pidx]["y"][0])
                    pointList += [pt]
                poly = geometry.Polygon([[p.x, p.y] for p in pointList])
                polyList += [poly]

            rects = prFramesAll[imgidx]["annorect"]
            prFramesAll[imgidx]["annorect"] = removeIgnoredPointsRects(rects,polyList)
            rects = gtFramesAll[imgidx]["annorect"]
            gtFramesAll[imgidx]["annorect"] = removeIgnoredPointsRects(rects,polyList)

    return gtFramesAll, prFramesAll

def generate_polygon(points):
  from shapely.geometry import Polygon, Point, LineString
  p = []
  p.append(points[0][1])
  for point in points:
    p.append(point[0])
  p.append(points[-1][1])
  for point in points[::-1]:
    p.append(point[2])
  eclipse_boundary = Polygon(p)

  p = []
  for point in points:
    p.append(point[1])
  center_line = LineString(p)

  return eclipse_boundary, center_line

def generate_polygon(points):
  from shapely.geometry import Polygon, Point, LineString
  p = []
  p.append(points[0][1])
  for point in points:
    p.append(point[0])
  p.append(points[-1][1])
  for point in points[::-1]:
    p.append(point[2])
  eclipse_boundary = Polygon(p)

  p = []
  for point in points:
    p.append(point[1])
  center_line = LineString(p)

  return eclipse_boundary, center_line

def draw(self, dc, f, **key):
        dc.SetPen(wx.Pen(Setting['color'], width=1, style=wx.SOLID))
        dc.SetTextForeground(Setting['tcolor'])
        font = wx.Font(10, wx.FONTFAMILY_DEFAULT, 
                       wx.FONTSTYLE_NORMAL, wx.FONTWEIGHT_NORMAL, False)

        dc.SetFont(font)
        dc.DrawLines([f(*i) for i in self.buf])
        for i in self.buf:dc.DrawCircle(f(*i),2)

        for pg in self.body:
            plg = Polygon(pg)
            dc.DrawLines([f(*i) for i in pg])
            for i in pg: dc.DrawCircle(f(*i),2)
            area, xy = plg.area, plg.centroid
            if self.unit!=None: 
                area *= self.unit[0]**2
            dc.DrawText('%.1f'%area, f(xy.x, xy.y))

def polygon_obb(polygon):
    '''
    Find the oriented bounding box of a Shapely polygon. 

    The OBB is always aligned with an edge of the convex hull of the polygon.

    Arguments
    -------------
    polygons: shapely.geometry.Polygon

    Returns
    -------------
    transform: (3,3) float, transformation matrix
               which will move input polygon from its original position 
               to the first quadrant where the AABB is the OBB
    extents:   (2,) float, extents of transformed polygon
    '''
    points = np.asanyarray(polygon.exterior.coords)
    return bounds.oriented_bounds_2D(points)

def polygon_hash(polygon):
    '''
    An approximate hash of a a shapely Polygon object.

    Arguments
    ---------
    polygon: shapely.geometry.Polygon object

    Returns
    ---------
    hash: (5) length list of hash representing input polygon
    '''
    result = [len(polygon.interiors),
              polygon.convex_hull.area,
              polygon.convex_hull.length,
              polygon.area, 
              polygon.length]
    return result

def random_polygon(segments=8, radius=1.0):
    '''
    Generate a random polygon with a maximum number of sides and approximate radius.

    Arguments
    ---------
    segments: int, the maximum number of sides the random polygon will have
    radius:   float, the approximate radius of the polygon desired

    Returns
    ---------
    polygon: shapely.geometry.Polygon object with random exterior, and no interiors. 
    '''
    angles = np.sort(np.cumsum(np.random.random(segments)*np.pi*2) % (np.pi*2))
    radii  = np.random.random(segments)*radius
    points = np.column_stack((np.cos(angles), np.sin(angles)))*radii.reshape((-1,1))
    points = np.vstack((points, points[0]))
    polygon = Polygon(points).buffer(0.0)
    if is_sequence(polygon):
        return polygon[0]
    return polygon

def clean_cut_polygon(polygon: Polygon) -> Polygon:
    interiors = []
    interiors.extend(cut_ring(polygon.exterior))
    exteriors = [(i, ring) for (i, ring) in enumerate(interiors) if ring.is_ccw]

    with suppress(AttributeError):
        delattr(polygon, 'c3nav_cache')

    if len(exteriors) != 1:
        raise ValueError('Invalid cut polygon!')
    exterior = interiors[exteriors[0][0]]
    interiors.pop(exteriors[0][0])

    for ring in polygon.interiors:
        interiors.extend(cut_ring(ring))

    return Polygon(exterior, interiors)

def _create_border(self, geometry: HybridGeometry, width, append=None):
        altitude = (np.vstack(chain(*(mesh.tolist() for mesh in geometry.faces)))[:, :, 2].max()+1)/1000
        geometry = self.buffered_bbox.intersection(geometry.geom)

        lines = tuple(chain(*(
            ((geom.exterior, *geom.interiors) if isinstance(geom, Polygon) else (geom,))
            for geom in getattr(geometry, 'geoms', (geometry,))
        )))

        if not lines:
            return np.empty((0, 3, 3+len(append)))

        lines = unary_union(lines).buffer(width, cap_style=CAP_STYLE.flat, join_style=JOIN_STYLE.mitre)

        vertices, faces = triangulate_polygon(lines)
        triangles = np.dstack((vertices[faces], np.full((faces.size, 1), fill_value=altitude).reshape((-1, 3, 1))))

        return self._append_to_vertices(triangles.astype(np.float32), append)

def _geometry_to_svg(self, geom):
        # scale and move geometry and create svg code for it
        if isinstance(geom, Polygon):
            return ('<path d="' +
                    ' '.join((('M %.1f %.1f L'+(' %.1f %.1f'*(len(ring.coords)-1))+' z') %
                              tuple((np.array(ring)*self.np_scale+self.np_offset).flatten()))
                             for ring in chain((geom.exterior,), geom.interiors))
                    + '"/>').replace('.0 ', ' ')
        if isinstance(geom, LineString):
            return (('<path d="M %.1f %.1f L'+(' %.1f %.1f'*(len(geom.coords)-1))+'"/>') %
                    tuple((np.array(geom)*self.np_scale+self.np_offset).flatten())).replace('.0 ', ' ')
        try:
            geoms = geom.geoms
        except AttributeError:
            return ''
        return ''.join(self._geometry_to_svg(g) for g in geoms)

def testLocator(self):
        image = mk_img(400, 600)

        # draw a rectangle
        A = (5, 5)
        B = (5, 300)
        C = (250, 5)
        D = (250, 300)
        ABCD = Polygon([A, B, D, C, A])
        image = draw_poly(image, ABCD)

        # locate it
        locator = BinaryLocator()
        located = locator.locate(image)
        polygons, labels = zip(*located)

        self.assertEqual(1, len(located), "One polygon found")
        self.assertTrue(ABCD.equals(polygons[0]), "Found polygon has the same shape")

        # test locate with an offset
        locator2 = BinaryLocator()
        located2 = locator2.locate(image, offset=(50, 40))
        polygons2, labels2 = zip(*located2)
        self.assertEqual(1, len(located2), "One polygon found")
        self.assertTrue(translate(ABCD, 50, 40).equals(polygons2[0]), "Found translated polygon")

def testLocator(self):
        image = mk_img(400, 600)

        # draw a rectangle
        A = (5, 80)
        B = (5, 300)
        C = (250, 80)
        D = (250, 300)
        ABCD = Polygon([A, B, D, C, A])
        image = draw_poly(image, ABCD)
        image, circle = draw_circle(image, 85, (500, 300), return_circle=True)

        # test locator
        locator = BinaryLocator()
        located = locator.locate(image)
        polygons, labels = zip(*located)

        self.assertEqual(2, len(polygons), "Two polygons found")
        self.assertTrue(ABCD.equals(polygons[1]), "Rectangle polygon is found")

        self.assertLessEqual(relative_error(polygons[0].area, np.pi * 85 * 85), 0.005)

def testLocate(self):
        image = mk_img(200, 300)

        # draw a rectangle
        A = (3, 40)
        B = (3, 150)
        C = (125, 40)
        D = (125, 150)
        ABCD = Polygon([A, B, D, C, A])
        image = draw_poly(image, ABCD, color=1)
        image = draw_circle(image, 40, (250, 150), color=2)

        # test locator
        locator = SemanticLocator(background=0)
        located = locator.locate(image)
        polygons, labels = zip(*located)

        self.assertEqual(2, len(polygons), "Two polygons found")
        self.assertTrue(ABCD.equals(polygons[0]), "Rectangle polygon is found")
        self.assertLessEqual(relative_error(polygons[1].area, np.pi * 40 * 40), 0.005)

def _get_intersecting_projects(search_polygon: Polygon, author_id: int):
        """ executes a database query to get the intersecting projects created by the author if provided """

        query = db.session.query(Project.id,
                                 Project.status,
                                 Project.default_locale,
                                 Project.geometry.ST_AsGeoJSON().label('geometry')) \
            .filter(ST_Intersects(Project.geometry,
                                  ST_MakeEnvelope(search_polygon.bounds[0],
                                                  search_polygon.bounds[1],
                                                  search_polygon.bounds[2],
                                                  search_polygon.bounds[3], 4326)))

        if author_id:
            query = query.filter(Project.author_id == author_id)

        return query.all()

def _geom_to_array(geom: BaseGeometry):
    if isinstance(geom, geometry.Point):
        yield np.array([(np.nan, GeomTypes.POINT)])
        yield np.asarray(geom.coords)
    elif isinstance(geom, geometry.LineString):
        yield np.array([(np.nan, GeomTypes.LINESTRING)])
        yield np.asarray(geom.coords)
    elif isinstance(geom, geometry.Polygon):
        for interior in geom.interiors:
            yield np.array([(np.nan, GeomTypes.POLYGON_HOLE)])
            yield np.asarray(interior)
        yield np.array([(np.nan, GeomTypes.POLYGON_SHELL)])
        yield np.asarray(geom.exterior)
    elif isinstance(geom, BaseMultipartGeometry):
        return chain.from_iterable(map(geom_to_array, geom))
    else:
        raise TypeError

def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        height = self.height
        width = self.width

        domain = Polygon([(0, 0), (0, height), (width, height), (width, 0)])
        target = LineString([(0, height / 2), (0, height)])
        spawn = Polygon([(0, 0),
                         (0, height / 2),
                         (width / 2, height / 2),
                         (width / 2, 0)])

        obstacles = LineString([(0, height / 2),
                                (width * self.ratio, height / 2)]) | \
                    domain.exterior - target

        self.obstacles = obstacles
        self.targets = [target]
        self.spawns = [spawn]
        self.domain = domain

def test_attributes():

    outl = ((0, 0), (0, 1), (1, 1.), (1, 0))
    r = Region_cls(1, 'Unit Square', 'USq', outl)  

    assert r.number == 1
    assert r.name == 'Unit Square'
    assert r.abbrev == 'USq'

    assert np.allclose(r.coords, outl)
    assert isinstance(r.coords, np.ndarray)

    assert r.polygon.equals(Polygon(outl))
    assert isinstance(r.polygon, Polygon)

    assert np.allclose(r.centroid, (0.5, 0.5))

def ellipses_intersect(ellipse1, ellipse2):
    a, b = ellipse_polyline((ellipse1, ellipse2))
    ea = LinearRing(a)
    eb = LinearRing(b)
    mp = ea.intersection(eb)
    #encloses = ea.contains(eb)
    #encloses = ea.contains(Point(ellipse2[0], ellipse2[1]))
    pa = Polygon(a)
    pb = Polygon(b)
    encloses = pa.contains(pb)

    x = [p.x for p in mp]
    y = [p.y for p in mp]
    if len(x) > 0:
        intersects = True
    else:
        intersects =  False

    return intersects, encloses

def get_way_centroid(self, way):
        """Calculate the centroid of a way

        way (object): overpy.Way object
        Returns dict of lat/lon
        """
        # Polygon has to have at least 3 points
        if len(way.nodes) >= 3:
            geom = []
            for node in way.nodes:
                geom.append((node.lon, node.lat))
            poly = Polygon(geom)
            cent = poly.centroid
            return {'lat': cent.y, 'lon': cent.x}
        elif len(way.nodes) == 2:
            # if way has 2 nodes, use average position
            lat = (way.nodes[0].lat + way.nodes[1].lat) / 2
            lon = (way.nodes[0].lon + way.nodes[1].lon) / 2
            return {'lat': lat, 'lon': lon}
        elif len(way.nodes) == 1:
            # if way has 1 node, use that position
            # (unusual and certainly a bug but possible)
            return {'lat': way.nodes[0].lat, 'lon': way.nodes[0].lon}
        else:
            raise RuntimeError

def lines_within_box(self, bounding_box):
        """
        Get all lines selected by a bounding box. Note that the whole line
        is selected when it partially falls into the box.
        :param bounding_box: the bounding coordinates in
            (minx, miny, maxx, maxy) or a Polygon instance
        :return: list, of LineInfo
        """
        subgraph = self.subgraph_within_box(bounding_box)
        lines = set()
        for edge in subgraph.edges():
            info = self.edge_info(edge)
            if info.line_index:
                lines.add(info.line_index)
        return [self.line_info(i) for i in lines]

def test_polygonal_min(self):
        polygon = Polygon([(0.0, 0.0), (0.0, 33.0), (33.0, 33.0), (33.0, 0.0), (0.0, 0.0)])
        result = self.tiled_rdd.polygonal_min(polygon, float)

        self.assertEqual(result, [0.0, 0.0])

def test_polygonal_max(self):
        polygon = Polygon([(1.0, 1.0), (1.0, 10.0), (10.0, 10.0), (10.0, 1.0)])
        result = self.tiled_rdd.polygonal_max(polygon, float)

        self.assertEqual(result, [1.0, 1.0])

def test_polygonal_sum(self):
        polygon = Polygon([(0.0, 0.0), (0.0, 33.0), (33.0, 33.0), (33.0, 0.0), (0.0, 0.0)])
        result = self.tiled_rdd.polygonal_sum(polygon, float)

        self.assertEqual(result, [96.0, 96.0])

def test_polygonal_mean(self):
        polygon = Polygon([(1.0, 1.0), (1.0, 10.0), (10.0, 10.0), (10.0, 1.0)])
        result = self.tiled_rdd.polygonal_mean(polygon)

        self.assertEqual(result, [1.0, 1.0])

def polygonal_mean(self, geometry):
        """Finds the mean of all of the values for each band that are contained within the given geometry.

        Args:
            geometry (shapely.geometry.Polygon or shapely.geometry.MultiPolygon or bytes): A
                Shapely ``Polygon`` or ``MultiPolygon`` that represents the area where the summary
                should be computed; or a WKB representation of the geometry.

        Returns:
            [float]
        """

        return self._process_polygonal_summary(geometry, self.srdd.polygonalMean)