Python shapely.geometry.MultiPoint() Examples

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 polygon_iou(list1, list2):
    """
    Intersection over union between two shapely polygons.
    """
    polygon_points1 = np.array(list1).reshape(4, 2)
    poly1 = Polygon(polygon_points1).convex_hull
    polygon_points2 = np.array(list2).reshape(4, 2)
    poly2 = Polygon(polygon_points2).convex_hull
    union_poly = np.concatenate((polygon_points1, polygon_points2))
    if not poly1.intersects(poly2):  # this test is fast and can accelerate calculation
        iou = 0
    else:
        try:
            inter_area = poly1.intersection(poly2).area
            # union_area = poly1.area + poly2.area - inter_area
            union_area = MultiPoint(union_poly).convex_hull.area
            if union_area == 0:
                return 0
            iou = float(inter_area) / union_area
        except shapely.geos.TopologicalError:
            print('shapely.geos.TopologicalError occured, iou set to 0')
            iou = 0
    return iou 

def shapely_to_paths(g):
    if isinstance(g, geometry.Point):
        return []
    elif isinstance(g, geometry.LineString):
        return [list(g.coords)]
    elif isinstance(g, (geometry.MultiPoint, geometry.MultiLineString, geometry.MultiPolygon, geometry.collection.GeometryCollection)):
        paths = []
        for x in g:
            paths.extend(shapely_to_paths(x))
        return paths
    elif isinstance(g, geometry.Polygon):
        paths = []
        paths.append(list(g.exterior.coords))
        for interior in g.interiors:
            paths.extend(shapely_to_paths(interior))
        return paths
    else:
        raise Exception('unhandled shapely geometry: %s' % type(g)) 

def shapely_to_paths(g):
    if isinstance(g, geometry.Point):
        return []
    elif isinstance(g, geometry.LineString):
        return [list(g.coords)]
    elif isinstance(g, (geometry.MultiPoint, geometry.MultiLineString, geometry.MultiPolygon, geometry.collection.GeometryCollection)):
        paths = []
        for x in g:
            paths.extend(shapely_to_paths(x))
        return paths
    elif isinstance(g, geometry.Polygon):
        paths = []
        paths.append(list(g.exterior.coords))
        for interior in g.interiors:
            paths.extend(shapely_to_paths(interior))
        return paths
    else:
        raise Exception('unhandled shapely geometry: %s' % type(g)) 

def _poly_area_approximation(way, nodes):
    """ Compute the approximated area of an irregular OSM polygon.
        see: https://arachnoid.com/area_irregular_polygon/
             https://gist.github.com/robinkraft/c6de2f988c9d3f01af3c
    """
    points = []
    for node in way['nd']:
        points.append([float(nodes[node['ref']]['lat']), float(nodes[node['ref']]['lon'])])

    approx = MultiPoint(points).convex_hull
    # http://openstreetmapdata.com/info/projections
    proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:3857'))

    converted_approximation = transform(proj, approx)

    return converted_approximation.area 

def get_geom_type(geom):
    """Returns the HoloViews geometry type.

    Args:
        geom: A shapely geometry

    Returns:
        A string representing type of the geometry.
    """
    from shapely.geometry import (
        Point, LineString, Polygon, Ring, MultiPoint, MultiPolygon, MultiLineString
    )
    if isinstance(geom, (Point, MultiPoint)):
        return 'Point'
    elif isinstance(geom, (LineString, MultiLineString)):
        return 'Line'
    elif isinstance(geom, Ring):
        return 'Ring'
    elif isinstance(geom, (Polygon, MultiPolygon)):
        return 'Polygon' 

def polygon_iou(list1, list2):
    """
    Intersection over union between two shapely polygons.
    """
    polygon_points1 = np.array(list1).reshape(4, 2)
    poly1 = Polygon(polygon_points1).convex_hull
    polygon_points2 = np.array(list2).reshape(4, 2)
    poly2 = Polygon(polygon_points2).convex_hull
    union_poly = np.concatenate((polygon_points1, polygon_points2))
    if not poly1.intersects(poly2):  # this test is fast and can accelerate calculation
        iou = 0
    else:
        try:
            inter_area = poly1.intersection(poly2).area
            # union_area = poly1.area + poly2.area - inter_area
            union_area = MultiPoint(union_poly).convex_hull.area
            if union_area == 0:
                return 1
            iou = float(inter_area) / union_area
        except shapely.geos.TopologicalError:
            print('shapely.geos.TopologicalError occured, iou set to 0')
            iou = 0
    return iou 

def polygon_iou(list1, list2):
    """
    Intersection over union between two shapely polygons.
    """
    polygon_points1 = np.array(list1).reshape(4, 2)
    poly1 = Polygon(polygon_points1).convex_hull
    polygon_points2 = np.array(list2).reshape(4, 2)
    poly2 = Polygon(polygon_points2).convex_hull
    union_poly = np.concatenate((polygon_points1, polygon_points2))
    if not poly1.intersects(poly2):  # this test is fast and can accelerate calculation
        iou = 0
    else:
        try:
            inter_area = poly1.intersection(poly2).area
            # union_area = poly1.area + poly2.area - inter_area
            union_area = MultiPoint(union_poly).convex_hull.area
            if union_area == 0:
                return 0
            iou = float(inter_area) / union_area
        except shapely.geos.TopologicalError:
            print('shapely.geos.TopologicalError occured, iou set to 0')
            iou = 0
    return iou 

def _round_multipoint_coords(mpt, precision):
    """
    Round the coordinates of a shapely MultiPoint to some decimal precision.

    Parameters
    ----------
    mpt : shapely.geometry.MultiPoint
        the MultiPoint to round the coordinates of
    precision : int
        decimal precision to round coordinates to

    Returns
    -------
    shapely.geometry.MultiPoint
    """
    return MultiPoint([_round_point_coords(pt, precision) for pt in mpt]) 

def test_cover_geometry_empty_geoms(tiler):
    """Empty geometries should return empty iterators."""
    assert not cover_geometry(tiler, geometry.Point(), 0) == True
    assert not cover_geometry(tiler, geometry.Point(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiPoint(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiPoint(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.LineString(), 0) == True
    assert not cover_geometry(tiler, geometry.LineString(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiLineString(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiLineString(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.Polygon(), 0) == True
    assert not cover_geometry(tiler, geometry.Polygon(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiPolygon(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiPolygon(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.GeometryCollection(), 0) == True
    assert not cover_geometry(tiler, geometry.GeometryCollection(), [0, 1]) == True 

def nearest_neighbor_within(others, point, max_distance):
    """Find nearest point among others up to a maximum distance.

    Args:
        others: a list of Points or a MultiPoint
        point: a Point
        max_distance: maximum distance to search for the nearest neighbor

    Returns:
        A shapely Point if one is within max_distance, None otherwise
    """
    search_region = point.buffer(max_distance)
    interesting_points = search_region.intersection(MultiPoint(others))

    if not interesting_points:
        closest_point = None
    elif isinstance(interesting_points, Point):
        closest_point = interesting_points
    else:
        distances = [point.distance(ip) for ip in interesting_points
                     if point.distance(ip) > 0]
        closest_point = interesting_points[distances.index(min(distances))]

    return closest_point 

def nearest_neighbor_within(others, point, max_distance):
    """Find nearest point among others up to a maximum distance.

    Args:
        others: a list of Points or a MultiPoint
        point: a Point
        max_distance: maximum distance to search for the nearest neighbor

    Returns:
        A shapely Point if one is within max_distance, None otherwise
    """
    search_region = point.buffer(max_distance)
    interesting_points = search_region.intersection(MultiPoint(others))

    if not interesting_points:
        closest_point = None
    elif isinstance(interesting_points, Point):
        closest_point = interesting_points
    else:
        distances = [point.distance(ip) for ip in interesting_points
                     if point.distance(ip) > 0]
        closest_point = interesting_points[distances.index(min(distances))]

    return closest_point 

def test_CheckTessellationInput(self):
        df = self.df_buildings
        df.loc[144, "geometry"] = Polygon([(0, 0), (0, 1), (1, 0)])
        df.loc[145, "geometry"] = MultiPoint([(0, 0), (1, 0)]).buffer(0.55)
        df.loc[146, "geometry"] = affinity.rotate(df.geometry.iloc[0], 12)
        check = mm.CheckTessellationInput(self.df_buildings)
        assert len(check.collapse) == 1
        assert len(check.split) == 1
        assert len(check.overlap) == 2

        check = mm.CheckTessellationInput(self.df_buildings, collapse=False)
        assert len(check.split) == 1
        assert len(check.overlap) == 2

        check = mm.CheckTessellationInput(self.df_buildings, split=False)
        assert len(check.collapse) == 1
        assert len(check.overlap) == 2

        check = mm.CheckTessellationInput(self.df_buildings, overlap=False)
        assert len(check.collapse) == 1
        assert len(check.split) == 1

        check = mm.CheckTessellationInput(self.df_buildings, shrink=0)
        assert len(check.collapse) == 0
        assert len(check.split) == 0
        assert len(check.overlap) == 4 

def to_geopandas(data, xdim, ydim, columns=[], geom='point'):
    """Converts list of dictionary format geometries to spatialpandas line geometries.

    Args:
        data: List of dictionaries representing individual geometries
        xdim: Name of x-coordinates column
        ydim: Name of y-coordinates column
        ring: Whether the data represents a closed ring

    Returns:
        A spatialpandas.GeoDataFrame version of the data
    """
    from geopandas import GeoDataFrame
    from shapely.geometry import (
        Point, LineString, Polygon, MultiPoint, MultiPolygon, MultiLineString
    )
    poly = any('holes' in d for d in data) or geom == 'Polygon'
    if poly:
        single_type, multi_type = Polygon, MultiPolygon
    elif geom == 'Line':
        single_type, multi_type = LineString, MultiLineString
    else:
        single_type, multi_type = Point, MultiPoint

    converted = defaultdict(list)
    for geom_dict in data:
        geom_dict = dict(geom_dict)
        geom = geom_from_dict(geom_dict, xdim, ydim, single_type, multi_type)
        for c, v in geom_dict.items():
            converted[c].append(v)
        converted['geometry'].append(geom)

    return GeoDataFrame(converted, columns=['geometry']+columns) 

def compute_intersections(lines, crs):
    import itertools
    inters = []
    for line1, line2 in itertools.combinations(lines, 2):
        if line1.intersects(line2):
            inter = line1.intersection(line2)
            if "Point" == inter.type:
                inters.append(inter)
            elif "MultiPoint" == inter.type:
                inters.extend([pt for pt in inter])
            elif "MultiLineString" == inter.type:
                multiLine = [line for line in inter]
                first_coords = multiLine[0].coords[0]
                last_coords = multiLine[len(multiLine) - 1].coords[1]
                inters.append(Point(first_coords[0], first_coords[1]))
                inters.append(Point(last_coords[0], last_coords[1]))
            elif "GeometryCollection" == inter.type:
                for geom in inter:
                    if "Point" == geom.type:
                        inters.append(geom)
                    elif "MultiPoint" == geom.type:
                        inters.extend([pt for pt in geom])
                    elif "MultiLineString" == geom.type:
                        multiLine = [line for line in geom]
                        first_coords = multiLine[0].coords[0]
                        last_coords = multiLine[len(multiLine) - 1].coords[1]
                        inters.append(Point(first_coords[0], first_coords[1]))
                        inters.append(Point(last_coords[0], last_coords[1]))

    geometry = inters
    df = gdf(geometry=geometry, crs=crs)
    return df 

def iloc(cls, dataset, index):
        from shapely.geometry import MultiPoint
        rows, cols = index

        data = dict(dataset.data)
        geom = data['geometry']

        if isinstance(geom, MultiPoint):
            if isscalar(rows) or isinstance(rows, slice):
                geom = geom[rows]
            elif isinstance(rows, (set, list)):
                geom = MultiPoint([geom[r] for r in rows])
        data['geometry'] = geom
        return data 

def generate_stop_points(chunks: List[LineString]) -> List[Point]:
    all_points = []

    # First point is the first node on the first line
    first_chunk = chunks[0].coords
    first_pt = [Point(p) for p in first_chunk][0]
    all_points.append(first_pt)

    # Then for all other points, we get from the end of each line
    for chunk in chunks:
        last_pt = [Point(p) for p in chunk.coords][-1]
        all_points.append(last_pt)

    # Now we need to convert to a Shapely object
    ap_ma = MultiPoint(all_points)

    # So we can reproject back out of equal area to 4326
    project = partial(
        pyproj.transform,
        pyproj.Proj(init='epsg:2163'),  # source coordinate system
        pyproj.Proj(init='epsg:4326'))  # destination coordinate system

    ap_ma_reproj = transform(project, ap_ma)  # apply projection

    # Final step will be to pull out all points back into a list
    return [p for p in ap_ma_reproj] 

def get_geom_point(x, sites):
    relevant_sites = sites.loc[zip(*[x.site_id, x.version])]
    return box(*MultiPoint(relevant_sites.geometry).bounds) 

def test_geometry_coords_rounding():
    # test the rounding of geometry coordinates
    precision = 3

    shape1 = Point(1.123456, 2.123456)
    shape2 = ox.utils_geo.round_geometry_coords(shape1, precision)

    shape1 = MultiPoint([(1.123456, 2.123456), (3.123456, 4.123456)])
    shape2 = ox.utils_geo.round_geometry_coords(shape1, precision)

    shape1 = LineString([(1.123456, 2.123456), (3.123456, 4.123456)])
    shape2 = ox.utils_geo.round_geometry_coords(shape1, precision)

    shape1 = MultiLineString(
        [
            [(1.123456, 2.123456), (3.123456, 4.123456)],
            [(11.123456, 12.123456), (13.123456, 14.123456)],
        ]
    )
    shape2 = ox.utils_geo.round_geometry_coords(shape1, precision)

    shape1 = Polygon([(1.123456, 2.123456), (3.123456, 4.123456), (6.123456, 5.123456)])
    shape2 = ox.utils_geo.round_geometry_coords(shape1, precision)

    shape1 = MultiPolygon(
        [
            Polygon([(1.123456, 2.123456), (3.123456, 4.123456), (6.123456, 5.123456)]),
            Polygon([(16.123456, 15.123456), (13.123456, 14.123456), (12.123456, 11.123456)]),
        ]
    )
    shape2 = ox.utils_geo.round_geometry_coords(shape1, precision) 

def round_geometry_coords(shape, precision):
    """
    Round the coordinates of a shapely geometry to some decimal precision.

    Parameters
    ----------
    shape : shapely.geometry.geometry
        the geometry to round the coordinates of; must be one of {Point,
        MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon}
    precision : int
        decimal precision to round coordinates to

    Returns
    -------
    shapely.geometry.geometry
    """
    if isinstance(shape, Point):
        return _round_point_coords(shape, precision)

    elif isinstance(shape, MultiPoint):
        return _round_multipoint_coords(shape, precision)

    elif isinstance(shape, LineString):
        return _round_linestring_coords(shape, precision)

    elif isinstance(shape, MultiLineString):
        return _round_multilinestring_coords(shape, precision)

    elif isinstance(shape, Polygon):
        return _round_polygon_coords(shape, precision)

    elif isinstance(shape, MultiPolygon):
        return _round_multipolygon_coords(shape, precision)

    else:
        raise TypeError(f"cannot round coordinates of unhandled geometry type: {type(shape)}") 

def multipoint():
    return geometry.MultiPoint([geometry.Point(0, 0), geometry.Point(1, 1)]) 

def create_ordered_line(feature_collection):
    ordered_collection = np.array(order_lines(feature_collection))
    return LineString(MultiPoint(ordered_collection)) 

def FindOrBuildDpa(dpas, options, grants):
  """Find or build DPA for simulation purpose.

  If several DPA, select the one with most grants around (but DPA simulation
  options always override the logic).
  """
  if options.dpa:
    dpa_kml_file = options.dpa_kml or None
    dpa = dpa_mgr.BuildDpa(options.dpa, None, portal_dpa_filename=dpa_kml_file)
  else:
    if not dpas:
      raise ValueError('Config file not defining a DPA and no --dpa option used.')
    if len(dpas) == 1:
      dpa = dpas[0]
    else:
      # Roughly find the DPA wth most CBSD within 200km
      all_cbsds = sgeo.MultiPoint([(g.longitude, g.latitude) for g in grants])
      num_cbsds_inside = [len(dpa.geometry.buffer(2.5).intersection(all_cbsds))
                          for dpa in dpas]
      dpa = dpas[np.argmax(num_cbsds_inside)]

  try: dpa.margin_db
  except AttributeError:
    dpa.margin_db = 'linear(1.5)'
  try: dpa.geometry
  except AttributeError:
    try: dpa_geometry = zones.GetCoastalDpaZones()[dpa.name]
    except KeyError: dpa_geometry = zones.GetPortalDpaZones(kml_path=dpa_kml_file)[dpa.name]
    dpa.geometry = dpa_geometry.geometry
  if options.dpa_builder:
    dpa.protected_points = dpa_builder.DpaProtectionPoints(
        dpa.name, dpa.geometry, options.dpa_builder)
  if options.margin_db:  # Override `movelistMargin` directive.
    try: dpa.margin_db = float(options.margin_db)
    except ValueError: dpa.margin_db = options.margin_db

  return dpa 

def CreateCbrsPlot(grants, dpa=None, tag='', color='g',
                   subplot=None, fig=None):
  """Plots the grants in a map along with other entities (DPA, ..).

  Args:
    grants: An iterable of |CbsdGrantInfo| grants.
    dpa: A |dpa_mgr.Dpa|  DPA.

  Returns:
    ax: the 'matplotlib.axes' object
  """
  if not fig:
    fig = plt.figure()
  if not subplot:
    subplot = 111
  ax = fig.add_subplot(subplot, projection=ccrs.PlateCarree())
  # Finds the bounding box of all geometries (DPA and CBSDs).
  box = sgeo.box(*sgeo.MultiPoint(
    [(grant.longitude, grant.latitude) for grant in grants]).bounds)
  if hasattr(dpa, 'geometry'):
    box = sgeo.box(*box.union(dpa.geometry).bounds)
  elif dpa is not None:
    box = sgeo.box(*box.union([(pt.longitude, pt.latitude)
                               for pt in dpa.protected_points]).bounds)
  box_margin = 0.1 # about 10km
  box = box.buffer(box_margin)
  # Plot geometries.
  ax.axis([box.bounds[0], box.bounds[2], box.bounds[1], box.bounds[3]])
  ax.coastlines()
  ax.stock_img()
  PlotGrants(ax, grants, color=color)
  if dpa is not None:
    PlotDpa(ax, dpa, color='m')
  ax.set_title('%sDPA: %s' % (tag, dpa.name))
  return ax, fig

#----------------------------------------------
# Useful statistical routines
# Estimate PDF form a set of random samples
#  - using a Smooth kernel estimator 

def test_degenerate_shapes(self):
    # Test all degenerate shapes (points, lines) have area zero
    points = sgeo.MultiPoint([(4, 59), (6, 59), (6, 61), (4, 61)])
    line = sgeo.LineString([(4, 59), (6, 59), (6, 61), (4, 61)])
    ring = sgeo.LinearRing([(4, 59), (6, 59), (6, 61), (4, 61)])
    self.assertEqual(utils.GeometryArea(points), 0)
    self.assertEqual(utils.GeometryArea(line), 0)
    self.assertEqual(utils.GeometryArea(ring), 0) 

def test_cover_geometry_multipoint3(tiler, mpt):
    """A MultiPoint geometry."""
    tiles = [tile for tile in cover_geometry(tiler, mpt, [3, 4])]
    assert len(tiles) == 1 

def test_cover_geometry_multipoint2(tiler, mpt):
    """A MultiPoint geometry."""
    tiles = [tile for tile in cover_geometry(tiler, mpt, 12)]
    assert len(tiles) == 9
    assert set(tiles) == {(973, 1203, 12), (974, 1203, 12), (975, 1203, 12),
                          (973, 1204, 12), (973, 1205, 12), (974, 1204, 12),
                          (975, 1204, 12), (974, 1205, 12), (975, 1205, 12)} 

def test_cover_geometry_multipoint1(tiler, mpt):
    """A MultiPoint geometry."""
    tiles = [tile for tile in cover_geometry(tiler, mpt, 4)]
    assert len(tiles) == 1 

def mpt():
    return geometry.MultiPoint([(-94.39453125, 15.908203125),
                                (-94.306640625, 15.908203125),
                                (-94.306640625, 15.8203125),
                                (-94.39453125, 15.8203125)]) 

def optimizealpha(points, max_iterations=10000):
    """
    Solve for the alpha parameter.

    Attempt to determine the alpha parameter that best wraps the given set of
    points in one polygon without dropping any points.

    Note:  If the solver fails to find a solution, a value of zero will be
    returned, which when used with the alphashape function will safely return a
    convex hull around the points.

    Args:

        points (list): an iterable container of points
        max_iterations (int): maximum number of iterations while finding the
            solution

    Returns:

        float: The optimized alpha parameter

    """
    # Convert to a shapely multipoint object if not one already
    if USE_GP and isinstance(points, geopandas.GeoDataFrame):
        points = points['geometry']
    if not isinstance(points, MultiPoint):
        points = MultiPoint(list(points))

    # Set the bounds
    lower = 0.

    # Ensure the upper limit bounds the solution
    upper = sys.float_info.max
    if _testalpha(points, upper):
        logging.error('the max float value does not bound the alpha '
                      'parameter solution')
        return 0.

    # Begin the bisection loop
    counter = 0
    while (upper - lower) > numpy.finfo(float).eps * 2:
        # Bisect the current bounds
        test_alpha = (upper + lower) * .5

        # Update the bounds to include the solution space
        if _testalpha(points, test_alpha):
            lower = test_alpha
        else:
            upper = test_alpha

        # Handle exceeding maximum allowed number of iterations
        counter += 1
        if counter > max_iterations:
            logging.warning('maximum allowed iterations reached while '
                            'optimizing the alpha parameter')
            break
    return lower