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

def merge_to_multi_polygon(feature_collection: str, dissolve: bool) -> geojson.MultiPolygon:
        """
        Merge all geometries to a single multipolygon
        :param feature_collection: geojson feature collection str containing features
        :param dissolve: flag for wther to to dissolve internal boundaries.
        :return: geojson.MultiPolygon
        """
        parsed_geojson = GridService._to_shapely_geometries(json.dumps(feature_collection))
        multi_polygon = GridService._convert_to_multipolygon(parsed_geojson)
        if dissolve:
            multi_polygon = GridService._dissolve(multi_polygon)
        aoi_multi_polygon_geojson =  geojson.loads(json.dumps(mapping(multi_polygon)))

        # validate the geometry
        if type(aoi_multi_polygon_geojson) is not geojson.MultiPolygon:
            raise InvalidGeoJson('Area Of Interest: geometry must be a MultiPolygon')

        is_valid_geojson = geojson.is_valid(aoi_multi_polygon_geojson)
        if is_valid_geojson['valid'] == 'no':
            raise InvalidGeoJson(f"Area of Interest: Invalid MultiPolygon - {is_valid_geojson['message']}")

        return aoi_multi_polygon_geojson

def _update_feature(clip_to_aoi: bool, feature: dict, new_shape) -> dict:
        """
        Updates the feature with the new shape, and splittable property
        :param clip_to_aoi: value for feature's splittable property
        :param feature: feature to be updated
        :param new_shape: new shape to use for feature
        :return:
        """
        if clip_to_aoi:
            # update the feature with the clipped shape
            if new_shape.geom_type == 'Polygon':
                # shapely may return a POLYGON rather than a MULTIPOLYGON if there is just one intersection area
                new_shape = MultiPolygon([new_shape])
            feature['geometry'] = mapping(new_shape)
            feature['properties']['x'] = None
            feature['properties']['y'] = None
            feature['properties']['zoom'] = None
            feature['properties']['splittable'] = False
        return feature

def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            image = GeoImage.__getitem__(self, geometry)
            return image
        else:
            result = DaskImage.__getitem__(self, geometry)
            image = super(GeoDaskWrapper, self.__class__).__new__(self.__class__,
                                                            result.dask, result.name, result.chunks,
                                                            result.dtype, result.shape)

            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            return image

def __new__(cls, op):
        assert isinstance(op, DaskMeta)
        self = super(IpeImage, cls).create(op)
        self._ipe_op = op
        if self.ipe.metadata["georef"] is None:
            tfm = RatPolyTransform.from_rpcs(self.ipe.metadata["rpcs"])
        else:
            tfm = AffineTransform.from_georef(self.ipe.metadata["georef"])
        img_md = self.ipe.metadata["image"]
        xshift = img_md["minTileX"]*img_md["tileXSize"]
        yshift = img_md["minTileY"]*img_md["tileYSize"]
        self.__geo_transform__ = tfm + (xshift, yshift)
        self.__geo_interface__ = mapping(self._reproject(wkt.loads(self.ipe.metadata["image"]["imageBoundsWGS84"])))
        minx = img_md["minX"] - xshift
        maxx = img_md["maxX"] - xshift
        miny = img_md["minY"] - yshift
        maxy = img_md["maxY"] - yshift

        return self[:, miny:maxy, minx:maxx]

def get_scan_coords():
        coordinates = mapping(conf.BOUNDARIES)['coordinates']
        coords = coordinates[0]
        markers = [{
                'type': 'scanarea',
                'coords': coords
            }]
        for blacklist in coordinates[1:]:
            markers.append({
                    'type': 'scanblacklist',
                    'coords': blacklist
                })
        return markers

def test_point(self):
        m = mapping(Point(0, 0))
        self.failUnlessEqual(m['type'], 'Point')
        self.failUnlessEqual(m['coordinates'], (0.0, 0.0))

def get_region_name(longitude, latitude):
    '''
        Function to get the region name given a point
    :param latitude: point's latitude
    :param longitude: point's longitude
    :return:
    '''
    point = 'POINT(' + str(longitude) + ' ' + str(latitude) + ')'
    query = session.query(Region).filter(func.ST_Contains(Region.regionboundary, point)).first()
    if query is None:
        return "Region not found"
    else:
        return {"name": query.name, "polygon": mapping(wkb.loads(bytes(query.regionboundary.data)))['coordinates']}

def getProjectedPatch(polygon, m, edgecolor, facecolor, lw=1., zorder=10):
    polyjson = mapping(polygon)
    tlist = []
    for sequence in polyjson['coordinates']:
        lon, lat = list(zip(*sequence))
        x, y = m(lon, lat)
        tlist.append(tuple(zip(x, y)))
    polyjson['coordinates'] = tuple(tlist)
    ppolygon = shape(polyjson)
    patch = PolygonPatch(ppolygon, facecolor=facecolor, edgecolor=edgecolor,
                         zorder=zorder, linewidth=lw, fill=True, visible=True)
    return patch

def to_geojson(self, instance=None) -> dict:
        result = OrderedDict((
            ('type', 'Feature'),
            ('properties', self.get_geojson_properties(instance=instance)),
            ('geometry', format_geojson(mapping(self.geometry), round=False)),
        ))
        original_geometry = getattr(self, 'original_geometry', None)
        if original_geometry:
            result['original_geometry'] = format_geojson(mapping(original_geometry), round=False)
        return result

def _serialize(self, geometry=True, simple_geometry=False, **kwargs):
        result = super()._serialize(simple_geometry=simple_geometry, **kwargs)
        if geometry:
            result['geometry'] = format_geojson(mapping(self.geometry), round=False)
        if simple_geometry:
            result['point'] = (self.level_id, ) + tuple(round(i, 2) for i in self.point.coords[0])
            if not isinstance(self.geometry, Point):
                minx, miny, maxx, maxy = self.geometry.bounds
                result['bounds'] = ((int(math.floor(minx)), int(math.floor(miny))),
                                    (int(math.ceil(maxx)), int(math.ceil(maxy))))
        return result

def details_display(self):
        result = super().details_display()
        result['geometry'] = format_geojson(mapping(self.geometry), round=False)
        return result

def _serialize(self, geometry=True, **kwargs):
        result = super()._serialize(geometry=geometry, **kwargs)
        result['width'] = float(str(self.width))
        result['height'] = float(str(self.height))
        if geometry:
            result['buffered_geometry'] = format_geojson(mapping(self.buffered_geometry))
        return result

def to_geojson(self, *args, **kwargs):
        result = super().to_geojson(*args, **kwargs)
        result['original_geometry'] = result['geometry']
        result['geometry'] = format_geojson(mapping(self.buffered_geometry))
        return result

def as_features(dct):
    for index, key in enumerate(dct):
        row = dct[key]
        gj = {
            'geometry': mapping(row['geom']),
            'properties': {
                'id': index,
                'from_label': key[0],
                'to_label': key[1],
                'measure': row['measure']},
        }
        yield gj

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 geometry_mapping(geom):
    from shapely.geometry import mapping

    if isinstance(geom, Iterable):
        return [mapping(geo) for geo in geom]
    elif isinstance(geom, Mapping):
        return {name: mapping(geo) for name, geo in geom.items()}
    else:
        return mapping(geom) if geom else geom

def write_segments_shp(segments_shp_path, road_projection, segments_with_data, schema):
    """Writes all segments to shapefile (both intersections and individual segments)
    :param segments_shp_path: Path to shapefile to write
    :param road_projection: Projection of road data
    :param segments_with_data: List of tuples containing segments and segment data
    :param schema: Schema to use for writing shapefile
    """
    with fiona.open(segments_shp_path, 'w', driver='ESRI Shapefile',
                    schema=schema, crs=road_projection) as output:
        for segment_with_data in segments_with_data:
            segment, data = segment_with_data
            output.write({
                'geometry': mapping(segment),
                'properties': data
            })

def __new__(cls, access_token=os.environ.get("DG_MAPS_API_TOKEN"),
                url="https://api.mapbox.com/v4/digitalglobe.nal0g75k/{z}/{x}/{y}.png",
                zoom=22, **kwargs):
        _tms_meta = TmsMeta(access_token=access_token, url=url, zoom=zoom, bounds=kwargs.get("bounds"))
        self = super(TmsImage, cls).create(_tms_meta)
        self._base_args = {"access_token": access_token, "url": url, "zoom": zoom}
        self._tms_meta = _tms_meta
        self.__geo_interface__ = mapping(box(*_tms_meta.bounds))
        self.__geo_transform__ = _tms_meta.__geo_transform__
        g = self._parse_geoms(**kwargs)
        if g is not None:
            return self[g]
        else:
            return self

def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            if self._tms_meta._bounds is None:
                return self.aoi(geojson=mapping(geometry), from_proj=self.proj)
            image = GeoImage.__getitem__(self, geometry)
            image._tms_meta = self._tms_meta
            return image
        else:
            result = super(TmsImage, self).__getitem__(geometry)
            image = super(TmsImage, self.__class__).__new__(self.__class__,
                                                            result.dask, result.name, result.chunks,
                                                            result.dtype, result.shape)

            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            image._tms_meta = self._tms_meta
            return image

def __getitem__(self, geometry):
        g = shape(geometry)

        bounds = ops.transform(self.__geo_transform__.rev, g).bounds
        try:
            assert g in self, "Image does not contain specified geometry {} not in {}".format(g.bounds, self.bounds)
        except AssertionError as ae:
            warnings.warn(ae.args)

        image = self._slice_padded(bounds)
        image.__geo_interface__ = mapping(g)
        return image

def boundary_polygon_from_file(filename):
    # TODO: This should be refactored and moved into datacube.utils.geometry
    import shapely.ops
    from shapely.geometry import shape, mapping
    with fiona.open(filename) as input_region:
        joined = shapely.ops.unary_union(list(shape(geom['geometry']) for geom in input_region))
        final = joined.convex_hull
        crs = CRS(input_region.crs_wkt)
        boundary_polygon = Geometry(mapping(final), crs)
    return boundary_polygon

def get_tile_geometry(path, origin_espg, tolerance=500):
    """ Calculate the data and tile geometry for sentinel-2 tiles """

    with rasterio.open(path) as src:

        # Get tile geometry
        b = src.bounds
        tile_shape = Polygon([(b[0], b[1]), (b[2], b[1]), (b[2], b[3]), (b[0], b[3]), (b[0], b[1])])
        tile_geojson = mapping(tile_shape)

        # read first band of the image
        image = src.read(1)

        # create a mask of zero values
        mask = image == 0.

        # generate shapes of the mask
        novalue_shape = shapes(image, mask=mask, transform=src.affine)

        # generate polygons using shapely
        novalue_shape = [Polygon(s['coordinates'][0]) for (s, v) in novalue_shape]

        if novalue_shape:

            # Make sure polygons are united
            # also simplify the resulting polygon
            union = cascaded_union(novalue_shape)

            # generates a geojson
            data_shape = tile_shape.difference(union)

            # If there are multipolygons, select the largest one
            if data_shape.geom_type == 'MultiPolygon':
                areas = {p.area: i for i, p in enumerate(data_shape)}
                largest = max(areas.keys())
                data_shape = data_shape[areas[largest]]

            # if the polygon has interior rings, remove them
            if list(data_shape.interiors):
                data_shape = Polygon(data_shape.exterior.coords)

            data_shape = data_shape.simplify(tolerance, preserve_topology=False)
            data_geojson = mapping(data_shape)

        else:
            data_geojson = tile_geojson

        # convert cooridnates to degrees
        return (to_latlon(tile_geojson, origin_espg), to_latlon(data_geojson, origin_espg))

def getNoSampleGrid(yespoints, xvar, yvar, dx, h1, h2):
    """Return the grid from which "no" pixels can successfully be sampled.

    Args:
        yespoints: Sequence of (x,y) points (meters) where
            landslide/liquefaction was observed.
        xvar: Numpy array of centers of columns of sampling grid.
        yvar: Numpy array of centers of rows of sampling grid.
        dx: Sampling resolution in x and y (meters).
        h1: Minimum buffer size for sampling non-hazard points.
        h2: Maximum buffer size for sampling non-hazard points.

    Returns:
        Grid of shape (len(yvar),len(xvar)) where 1's represent pixels from
        which "no" values can be sampled.
    """

    shp = (len(xvar), len(yvar))
    west = xvar.min() - dx / 2.0  # ??
    north = yvar.max() + dx / 2.0  # ??
    affine = affine_from_corner(west, north, dx, dx)
    donuts = []
    holes = []
    for h, k in yespoints:
        donut = createCirclePolygon(h, k, h2, dx)
        hole = createCirclePolygon(h, k, h1, dx)
        donuts.append(donut)
        holes.append(hole)
    donutburn = ((mapping(g), 1) for g in donuts)
    holeburn = ((mapping(g), 2) for g in holes)
    # we only want those pixels set where the polygon encloses the center point
    alltouched = False
    donutimg = rasterio.features.rasterize(donutburn,
                                           out_shape=shp,
                                           transform=affine,
                                           all_touched=alltouched)
    holeimg = rasterio.features.rasterize(holeburn,
                                          out_shape=shp,
                                          transform=affine,
                                          all_touched=alltouched)
    holeimg[holeimg == 0] = 1
    holeimg[holeimg == 2] = 0
    sampleimg = np.bitwise_and(donutimg, holeimg)
    return sampleimg

def getPolygons(self, keepUncompletePolygons=True):
        features = []

        # for polygons in relations (especially outer/innner) 
        # online forums suggest to not trust the outer tag
        # so we rather test it by ourselves
        for k, members in self.relations.iteritems():
            polygons = []
            properties = {}
            for osmid,relType,role in members:
                # if the way does not belong to a building
                if osmid not in self.ways:
                    continue
                elif not self.ways[osmid][1]:
                    continue
                polygon = Polygon([self.coordinates[ref] for ref in self.ways[osmid][0]])
                polygons.append(polygon)

            if len(polygons) < 2:
                continue

            # sort polygons in whithinness order
            # if they are all inside the first we make a polygon with holes
            polygons = sorted(polygons, key=Within, reverse=True)   
            if Within(polygons[1]) < Within(polygons[0]):
                polygon = geojson.Polygon([list(mapping(p)['coordinates'][0]) for p in polygons])
                refs,building,amenity,wheelchair,buildingLevels = self.ways[members[0][0]]
                properties = {'osmid':k,'amenity':amenity,'wheelchair':wheelchair,'levels':buildingLevels}
                features.append(geojson.Feature(geometry=polygon, properties=properties))

                for osmid,relType,role in members:
                    del self.ways[osmid]

        for osmid, (refs,building,amenity,wheelchair,buildingLevels) in self.ways.iteritems():
            if building:
                polygon = [[tuple(self.coordinates[ref]) for ref in refs if ref in self.coordinates]] 
                if keepUncompletePolygons:
                    if len(polygon[0]) < 3: 
                        continue
                elif len(polygon[0]) != len(refs):
                    continue
                properties = {'osmid':osmid,'amenity':amenity,'wheelchair':wheelchair,'levels':buildingLevels}
                features.append(geojson.Feature(geometry=geojson.Polygon(polygon), properties=properties))
        featureCollection = geojson.FeatureCollection(features)
        return featureCollection

def test_getPolygons(self, keepUncompletePolygons=True):
        features = []

        # for polygons in relations (especially outer/innner) 
        # online forums suggest to not trust the outer tag
        # so we rather test it by ourselves
        for k, members in self.relations.iteritems():
            polygons = []
            properties = {}
            for osmid,relType,role in members:
                # if the way does not belong to a building
                if osmid not in self.ways:
                    continue
                elif not self.ways[osmid][1]:
                    continue
                polygon = Polygon([self.coordinates[ref] for ref in self.ways[osmid][0]])
                polygons.append(polygon)

            if len(polygons) < 2:
                continue

            # sort polygons in whithinness order
            # if they are all inside the first we make a polygon with holes
            polygons = sorted(polygons, key=Within, reverse=True)   
            if Within(polygons[1]) < Within(polygons[0]):
                polygon = geojson.Polygon([list(mapping(p)['coordinates'][0]) for p in polygons])
                refs,building,amenity,wheelchair,buildingLevels = self.ways[members[0][0]]
                properties = {'osmid':k,'amenity':amenity,'wheelchair':wheelchair,'levels':buildingLevels}
                features.append(geojson.Feature(geometry=polygon, properties=properties))

                for osmid,relType,role in members:
                    del self.ways[osmid]

        for osmid, (refs,building,amenity,wheelchair,buildingLevels) in self.ways.iteritems():
            if building:
                polygon = [[tuple(self.coordinates[ref]) for ref in refs if ref in self.coordinates]] 
                if keepUncompletePolygons:
                    if len(polygon[0]) < 3: 
                        continue
                elif len(polygon[0]) != len(refs):
                    continue
                properties = {'osmid':osmid,'amenity':amenity,'wheelchair':wheelchair,'levels':buildingLevels}
                features.append(geojson.Feature(geometry=geojson.Polygon(polygon), properties=properties))
        featureCollection = geojson.FeatureCollection(features)
        return featureCollection