python类MultiPolygon()的实例源码

def parse_poly(lines):
    """ Parse an Osmosis/Google polygon filter file.

    :Params:
        lines : open file pointer

    :Returns:
        shapely.geometry.MultiPolygon object.

    .. note::
        http://wiki.openstreetmap.org/wiki/Osmosis/Polygon_Filter_File_Format

    """
    in_ring = False
    coords = []

    for index, line in enumerate(lines):
        if index == 0:
            # ignore meta/garbage.
            continue
        elif index == 1:
            coords.append([[], []])
            ring = coords[-1][0]
            in_ring = True
        elif in_ring and line.strip() == 'END':
            in_ring = False
        elif in_ring:
            ring.append(map(float, line.split()))
        elif not in_ring and line.strip() == 'END':
            break
        elif not in_ring and line.startswith('!'):
            coords[-1][1].append([])
            ring = coords[-1][1][-1]
            in_ring = True
        elif not in_ring:
            coords.append([[], []])
            ring = coords[-1][0]
            in_ring = True

    return MultiPolygon(coords)

def parse_poly(lines):
    """ Parse an Osmosis/Google polygon filter file.

    :Params:
        lines : open file pointer

    :Returns:
        shapely.geometry.MultiPolygon object.

    .. note::
        http://wiki.openstreetmap.org/wiki/Osmosis/Polygon_Filter_File_Format

    """
    in_ring = False
    coords = []

    for index, line in enumerate(lines):
        if index == 0:
            # ignore meta/garbage.
            continue
        elif index == 1:
            coords.append([[], []])
            ring = coords[-1][0]
            in_ring = True
        elif in_ring and line.strip() == 'END':
            in_ring = False
        elif in_ring:
            ring.append(map(float, line.split()))
        elif not in_ring and line.strip() == 'END':
            break
        elif not in_ring and line.startswith('!'):
            coords[-1][1].append([])
            ring = coords[-1][1][-1]
            in_ring = True
        elif not in_ring:
            coords.append([[], []])
            ring = coords[-1][0]
            in_ring = True

    return MultiPolygon(coords)

def classes(self):
        return {
            'polygon': (Polygon, ),
            'multipolygon': (Polygon, MultiPolygon),
            'linestring': (LineString, ),
            'point': (Point, )
        }[self.geomtype]

def triangulate_polygon(geometry: Union[Polygon, MultiPolygon], keep_holes=False):
    if isinstance(geometry, Polygon):
        return _triangulate_polygon(geometry, keep_holes=keep_holes)

    vertices = deque()
    faces = deque()

    offset = 0
    for polygon in geometry.geoms:
        new_vertices, new_faces = _triangulate_polygon(polygon, keep_holes=keep_holes)
        vertices.append(new_vertices)
        faces.append(new_faces+offset if offset else new_faces)
        offset += len(new_vertices)

    return np.vstack(vertices), np.vstack(faces)

def clean_geometry(geometry):
    """
    if the given geometry is a Polygon and invalid, try to make it valid if it results in a Polygon (not MultiPolygon)
    """
    if geometry.is_valid:
        return geometry

    if isinstance(geometry, Polygon):
        return geometry.buffer(0)

    return geometry

def assert_multipolygon(geometry: Union[Polygon, MultiPolygon, GeometryCollection]) -> List[Polygon]:
    """
    given a Polygon or a MultiPolygon, return a list of Polygons
    :param geometry: a Polygon or a MultiPolygon
    :return: a list of Polygons
    """
    if geometry.is_empty:
        return []
    if isinstance(geometry, Polygon):
        return [geometry]
    return [geom for geom in geometry.geoms if isinstance(geom, Polygon)]

def shapely_to_mpl(geometry):
    """
    convert a shapely Polygon or Multipolygon to a matplotlib Path
    :param polygon: shapely Polygon or Multipolygon
    :return: MplPathProxy
    """
    if isinstance(geometry, Polygon):
        return MplPolygonPath(geometry)
    elif isinstance(geometry, MultiPolygon) or geometry.is_empty or isinstance(geometry, GeometryCollection):
        return MplMultipolygonPath(geometry)
    raise TypeError

def _create_square(x, y, zoom) -> geojson.MultiPolygon:
        """
        Function for creating a geojson.MultiPolygon square representing a single OSM tile grid square
        :param x: osm tile grid x
        :param y: osm tile grid y
        :param zoom: osm tile grid zoom level
        :return: geojson.MultiPolygon in EPSG:4326
        """
        # Maximum resolution
        MAXRESOLUTION = 156543.0339

        # X/Y axis limit
        max = MAXRESOLUTION * 256 / 2

        # calculate extents
        step = max / (2 ** (zoom - 1))
        xmin = x * step - max
        ymin = y * step - max
        xmax = (x + 1) * step - max
        ymax = (y + 1) * step - max

        # make a shapely multipolygon
        multipolygon = MultiPolygon([Polygon([(xmin, ymin), (xmax, ymin),
                                              (xmax, ymax), (xmin, ymax)])])

        # use the database to transform the geometry from 3857 to 4326
        transformed_geometry = ST_Transform(shape.from_shape(multipolygon, 3857), 4326)

        # use DB to get the geometry as geojson
        return geojson.loads(db.engine.execute(transformed_geometry.ST_AsGeoJSON()).scalar())

def trim_grid_to_aoi(grid_dto: GridDTO) -> geojson.FeatureCollection:
        """
        Removes grid squares not intersecting with the aoi.  Optionally leaves partialy intersecting task squares
        complete or clips them exactly to the AOI outline
        :param grid_dto: the dto containing
        :return: geojson.FeatureCollection trimmed task grid
        """
        # get items out of the dto
        grid = grid_dto.grid
        aoi = grid_dto.area_of_interest
        clip_to_aoi = grid_dto.clip_to_aoi

        # create a shapely shape from the aoi
        aoi_multi_polygon_geojson = GridService.merge_to_multi_polygon(aoi, dissolve=True)

        aoi_multi_polygon = shapely.geometry.shape(aoi_multi_polygon_geojson)
        intersecting_features = []
        for feature in grid['features']:
            # create a shapely shape for the tile
            tile = shapely.geometry.shape(feature['geometry'])
            if aoi_multi_polygon.contains(tile):
                # tile is completely within aoi, use as is
                intersecting_features.append(feature)
            else:
                intersection = aoi_multi_polygon.intersection(tile)
                if intersection.is_empty or intersection.geom_type not in ['Polygon', 'MultiPolygon']:
                    continue  # this intersections which are not polygons or which are completely outside aoi
                # tile is partially intersecting the aoi
                clipped_feature = GridService._update_feature(clip_to_aoi, feature, intersection)
                intersecting_features.append(clipped_feature)
        return geojson.FeatureCollection(intersecting_features)

def tasks_from_aoi_features(feature_collection: str) -> geojson.FeatureCollection:
        """
        Creates a geojson feature collection of tasks from an aoi feature collection
        :param feature_collection:
        :return: task features
        """
        parsed_geojson = GridService._to_shapely_geometries(json.dumps(feature_collection))
        tasks = []
        for feature in parsed_geojson:
            if not isinstance(feature.geometry, MultiPolygon):
                feature.geometry = MultiPolygon([feature.geometry])
            # put the geometry back to geojson

            if feature.geometry.has_z:
                # Strip Z dimension, as can't persist geometry otherewise.  Most likely exists in KML data
                feature.geometry = shapely.ops.transform(GridService._to_2d, feature.geometry)

            feature.geometry = shapely.geometry.mapping(feature.geometry)

            # set default properties
            feature.properties = {
                'x': None,
                'y': None,
                'zoom': None,
                'splittable': False
            }

            tasks.append(feature)

        return geojson.FeatureCollection(tasks)

def _dissolve(geoms: MultiPolygon) -> MultiPolygon:
        """
        dissolves a Multipolygons
        :return: Multipolygon
        """
        # http://toblerity.org/shapely/manual.html#shapely.ops.cascaded_union
        geometry = cascaded_union(geoms)
        if geometry.geom_type == 'Polygon':
            # shapely may return a POLYGON rather than a MULTIPOLYGON if there is just one shape
            # force Multipolygon
            geometry = MultiPolygon([geometry])
        return geometry

def polygons(self):
        """list of shapely Polygon/ MultiPolygon of the regions"""
        return self.combiner('polygon')

def _is_polygon(self):
        """is there at least one region that was a Polygon/ MultiPolygon

        ."""
        return np.any(np.array(self.combiner('_is_polygon')))


# add the plotting method

def polygon(self):
        """shapely Polygon or MultiPolygon of the region"""

        if self._polygon is None:
            self._polygon = Polygon(self.coords)
        return self._polygon

def test_polygon_input():

    # polygon closes open paths
    outl1 = ((0, 0), (0, 1), (1, 1.), (1, 0), (0, 0))
    outl2 = ((1, 1), (1, 2), (2, 2.), (2, 1), (1, 1))
    nan = np.ones(shape=(1, 2)) * np.nan
    outl = np.vstack((outl1, nan, outl2))

    outl_poly = MultiPolygon([Polygon(outl1), Polygon(outl2)])

    r = Region_cls(1, 'Unit Square', 'USq', outl_poly)

    assert np.allclose(r.coords, outl, equal_nan=True)
    assert r.polygon == outl_poly

def features_on_raster(self, filtercodes=None, raster_extent=None):
        from shapely.geometry import MultiPolygon
        raster_extent = raster_extent if raster_extent else self.demDataExtent

        if filtercodes:
            features = [feat for feat in self.features if self.fieldName_Description in feat['fields'] and feat['fields'][self.fieldName_Description] in filtercodes]
        else:
            features = self.features
        for feat in features:
            if not MultiPolygon(raster_extent).contains(feat['geometry']):
                return False
        return True

def singlepart_features_test(self):
        if self.features:
            if any(feat['geometry'].type in ['MultiLineString', 'MultiPolygon', 'MultiPoint'] for feat in self.features):
                return False
        return True

def start_stop_on_raster(self, raster_extent=None):
        from shapely.geometry import Point, MultiPolygon
        raster_extent = raster_extent if raster_extent else self.demDataExtent
        if self.features and raster_extent:
            for feat in self.features:
                if not (MultiPolygon(raster_extent).contains(Point(feat['geometry'].coords[0])) and
                        MultiPolygon(raster_extent).contains(Point(feat['geometry'].coords[-1]))):
                    return False
            return True

def _process_polygonal_summary(geometry, operation):
        if isinstance(geometry, (Polygon, MultiPolygon)):
            geometry = shapely.wkb.dumps(geometry)
        if not isinstance(geometry, bytes):
            raise TypeError("Expected geometry to be bytes but given this instead", type(geometry))

        return list(operation(geometry))

def mask_to_poly(image_id):
    global size
    xx = []
    for version, epoch in nets:
        x = joblib.load('/data/raw_preds/{}-{}/{}.pkl'.format(version, epoch, image_id))
        x = cv2.resize(x.transpose((1, 2, 0)), (size, size))
        xx.append(x)

    xx = np.stack(xx, axis=3)
    preds = np.mean(xx, axis=3)
#    preds = joblib.load('raw_preds/{}-{}/{}.pkl'.format('v63', 2700, image_id))
#    preds = cv2.resize(preds.transpose((1, 2, 0)), (size, size))
    preds = preds.transpose((2, 0, 1)).copy()
    print(image_id)

    if n_out == 10:
#        preds = (preds > 0.3).astype(np.uint8)

        thresholds = np.array([0.4, 0.3, 0.3, 0.3, 0.7,
                               0.4, 0.4, 0.4, 0.04, 0.04]).reshape((10, 1))
        preds = (preds.reshape((10, -1)) > thresholds).reshape((10, size, size))
        preds = preds.astype(np.uint8)
    else:
        preds = np.argmax(preds, axis=0)
        preds = unsoft(preds)

    rg = colorize_raster(preds.transpose((1, 2, 0)))
    rg_size = 700
    rg = cv2.resize(rg, (rg_size, rg_size))
    im = get_rgb_image(image_id, rg_size, rg_size)
    rg = np.hstack([rg, im])
    cv2.imwrite('raw_blend_temp5/{}.png'.format(image_id), rg)
    joblib.dump(preds.astype(np.float32), 'raw_preds/raw_blend5/{}.pkl'.format(image_id))
    return

    shs = []
    for i in range(10):
        mask = preds[i]

        y_sf, x_sf = get_scale_factor(image_id, mask.shape[0], mask.shape[1])
        y_sf = 1. / y_sf
        x_sf = 1. / x_sf

        sh = polygonize_cv(mask)
#        sh = polygonize_sk((mask>0)*255, 0)
#        sh = (sh1.buffer(0).intersection(sh2.buffer(0))).buffer(0)

#        if not sh.is_valid:
#            sh = sh.buffer(0)
        sh = affinity.scale(sh, xfact=x_sf, yfact=y_sf, origin=(0, 0, 0))

        try:
            sh = MultiPolygon(sh)
        except:
            print('ERRRROR!!')
            sh = MultiPolygon()

        shs.append(sh)

    return shs