Python shapely.ops方法代码示例

# 需要导入模块: import shapely [as 别名]
# 或者: from shapely import ops [as 别名]
def mask_to_poly(mask, min_polygon_area_th=MIN_AREA):
    shapes = rasterio.features.shapes(mask.astype(np.int16), mask > 0)
    poly_list = []
    mp = shapely.ops.cascaded_union(
        shapely.geometry.MultiPolygon([
            shapely.geometry.shape(shape)
            for shape, value in shapes
        ]))

    if isinstance(mp, shapely.geometry.Polygon):
        df = pd.DataFrame({
            'area_size': [mp.area],
            'poly': [mp],
        })
    else:
        df = pd.DataFrame({
            'area_size': [p.area for p in mp],
            'poly': [p for p in mp],
        })

    df = df[df.area_size > min_polygon_area_th].sort_values(
        by='area_size', ascending=False)
    df.loc[:, 'wkt'] = df.poly.apply(lambda x: shapely.wkt.dumps(
        x, rounding_precision=0))
    df.loc[:, 'bid'] = list(range(1, len(df) + 1))
    df.loc[:, 'area_ratio'] = df.area_size / df.area_size.max()
    return df 

# 需要导入模块: import shapely [as 别名]
# 或者: from shapely import ops [as 别名]
def iter_data(self, geom_type, field_indices, slicing, mask_poly, mask_rect, clip):
        clip_poly = None
        if mask_poly is not None:
            mask_poly = conv.ogr_of_shapely(mask_poly)
            if clip:
                clip_poly = mask_poly
        elif mask_rect is not None:
            if clip:
                clip_poly = conv.ogr_of_shapely(sg.box(*mask_rect))

        ftr = None # Necessary to prevent the old swig bug
        geom = None # Necessary to prevent the old swig bug
        with self.acquire_driver_object() as (_, lyr):
            for ftr in self.iter_features_driver(slicing, mask_poly, mask_rect, lyr):
                geom = ftr.geometry()
                if geom is None or geom.IsEmpty():
                    # `geom is None` and `geom.IsEmpty()` is not exactly the same case, but whatever
                    geom = None
                    if not self.back_ds.allow_none_geometry: # pragma: no cover
                        raise Exception(
                            'None geometry in feature '
                            '(Set `allow_none_geometry=True` in Dataset constructor to silence)'
                        )
                else:
                    if clip:
                        geom = geom.Intersection(clip_poly)
                        assert not geom.IsEmpty()
                    geom = conv.shapely_of_ogr(geom)
                    if self.to_work:
                        geom = shapely.ops.transform(self.to_work, geom)
                    if geom_type == 'coordinates':
                        geom = sg.mapping(geom)['coordinates']
                    elif geom_type == 'geojson':
                        geom = sg.mapping(geom)

                yield (geom,) + tuple(
                    self._type_of_field_index[index](ftr.GetField(index))
                    if ftr.GetField(index) is not None
                    else None
                    for index in field_indices
                )

        # Necessary to prevent the old swig bug
        # https://trac.osgeo.org/gdal/ticket/6749
        del geom
        del ftr
        del clip_poly
        del mask_rect, mask_poly 

# 需要导入模块: import shapely [as 别名]
# 或者: from shapely import ops [as 别名]
def _ogr_of_geom(self, geom, geom_type):
        if geom_type is None:
            geom = geom
        elif self.to_virtual:
            if geom_type == 'coordinates':
                geom = sg.shape({
                    'type': self.type,
                    'coordinates': geom,
                })
            geom = shapely.ops.transform(self.to_virtual, geom)
            # geom = conv.ogr_of_shapely(geom)
            # TODO: Use json and unit test
            mapping = sg.mapping(geom)
            geom = conv.ogr_of_coordinates(
                mapping['coordinates'],
                mapping['type'],
            )
            if geom is None: # pragma: no cover
                raise ValueError('Could not convert `{}` of type `{}` to `ogr.Geometry`'.format(
                    geom_type, self.type
                ))
        elif geom_type == 'coordinates':
            geom = conv.ogr_of_coordinates(geom, self.type)
            if geom is None: # pragma: no cover
                raise ValueError('Could not convert `{}` of type `{}` to `ogr.Geometry`'.format(
                    geom_type, self.type
                ))
        elif geom_type == 'shapely':
            # geom = conv.ogr_of_shapely(geom)
            # TODO: Use json and unit test
            mapping = sg.mapping(geom)
            geom = conv.ogr_of_coordinates(
                mapping['coordinates'],
                mapping['type'],
            )
            if geom is None: # pragma: no cover
                raise ValueError('Could not convert `{}` of type `{}` to `ogr.Geometry`'.format(
                    geom_type, self.type
                ))
        else:
            assert False # pragma: no cover
        return geom

    # Misc ************************************************************************************** ** 

# 需要导入模块: import shapely [as 别名]
# 或者: from shapely import ops [as 别名]
def _create_merged_climate_zones_metadata(county_metadata):
    from shapely.ops import cascaded_union

    iecc_climate_zone_polygons = defaultdict(list)
    iecc_moisture_regime_polygons = defaultdict(list)
    ba_climate_zone_polygons = defaultdict(list)

    for county in county_metadata.values():
        polygon = county["polygon"]
        iecc_climate_zone = county.get("iecc_climate_zone")
        iecc_moisture_regime = county.get("iecc_moisture_regime")
        ba_climate_zone = county.get("ba_climate_zone")
        if iecc_climate_zone is not None:
            iecc_climate_zone_polygons[iecc_climate_zone].append(polygon)
        if iecc_moisture_regime is not None:
            iecc_moisture_regime_polygons[iecc_moisture_regime].append(polygon)
        if ba_climate_zone is not None:
            ba_climate_zone_polygons[ba_climate_zone].append(polygon)

    iecc_climate_zone_metadata = {}
    for iecc_climate_zone, polygons in iecc_climate_zone_polygons.items():
        polygon = cascaded_union(polygons)
        polygon = polygon.simplify(0.01)
        iecc_climate_zone_metadata[iecc_climate_zone] = {
            "iecc_climate_zone": iecc_climate_zone,
            "polygon": polygon,
            "geometry": to_geojson(polygon),
        }

    iecc_moisture_regime_metadata = {}
    for iecc_moisture_regime, polygons in iecc_moisture_regime_polygons.items():
        polygon = cascaded_union(polygons)
        polygon = polygon.simplify(0.01)
        iecc_moisture_regime_metadata[iecc_moisture_regime] = {
            "iecc_moisture_regime": iecc_moisture_regime,
            "polygon": polygon,
            "geometry": to_geojson(polygon),
        }

    ba_climate_zone_metadata = {}
    for ba_climate_zone, polygons in ba_climate_zone_polygons.items():
        polygon = cascaded_union(polygons)
        polygon = polygon.simplify(0.01)
        ba_climate_zone_metadata[ba_climate_zone] = {
            "ba_climate_zone": ba_climate_zone,
            "polygon": polygon,
            "geometry": to_geojson(polygon),
        }

    return (
        iecc_climate_zone_metadata,
        iecc_moisture_regime_metadata,
        ba_climate_zone_metadata,
    ) 

# 需要导入模块: import shapely [as 别名]
# 或者: from shapely import ops [as 别名]
def alpha_shape(self, coords, alpha):
        import shapely.geometry as geometry
        from shapely.ops import cascaded_union, polygonize
        from scipy.spatial import Delaunay
        """
        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] ])
        
        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