Python rasterio.features.shapes() Examples

def union(inputtiles, parsenames):

    tiles = sutils.tile_parser(inputtiles, parsenames)

    xmin, xmax, ymin, ymax = sutils.get_range(tiles)

    zoom = sutils.get_zoom(tiles)

    # make an array of shape (xrange + 3, yrange + 3)
    burn = sutils.burnXYZs(tiles, xmin, xmax, ymin, ymax, 0)

    nw = mercantile.xy(*mercantile.ul(xmin, ymin, zoom))

    se = mercantile.xy(*mercantile.ul(xmax + 1, ymax + 1, zoom))

    aff = Affine(((se[0] - nw[0]) / float(xmax - xmin + 1)), 0.0, nw[0],
        0.0, -((nw[1] - se[1]) / float(ymax - ymin + 1)), nw[1])

    unprojecter = sutils.Unprojecter()

    unionedTiles = [
        {
            'geometry': unprojecter.unproject(feature),
            'properties': {},
            'type': 'Feature'
        } for feature, shapes in features.shapes(np.asarray(np.flipud(np.rot90(burn)).astype(np.uint8), order='C'), transform=aff) if shapes == 1
    ]

    return unionedTiles 

def crowns_to_polys_raster(self):
        ''' Converts tree crown raster to individual polygons and stores them
        in the tree dataframe
        '''
        polys = []
        for feature in rioshapes(self.crowns, mask=self.crowns.astype(bool)):

            # Convert pixel coordinates to lon/lat
            edges = feature[0]['coordinates'][0].copy()
            for i in range(len(edges)):
                edges[i] = self._to_lonlat(*edges[i], self.resolution)

            # poly_smooth = self.smooth_poly(Polygon(edges), s=None, k=9)
            polys.append(Polygon(edges))
        self.trees.crown_poly_raster = polys 

def _compute_image_stats_chunked(dataset: 'DatasetReader') -> Optional[Dict[str, Any]]:
        """Compute statistics for the given rasterio dataset by looping over chunks."""
        from rasterio import features, warp, windows
        from shapely import geometry

        total_count = valid_data_count = 0
        tdigest = TDigest()
        sstats = SummaryStats()
        convex_hull = geometry.Polygon()

        block_windows = [w for _, w in dataset.block_windows(1)]

        for w in block_windows:
            with warnings.catch_warnings():
                warnings.filterwarnings('ignore', message='invalid value encountered.*')
                block_data = dataset.read(1, window=w, masked=True)

            total_count += int(block_data.size)
            valid_data = block_data.compressed()

            if valid_data.size == 0:
                continue

            valid_data_count += int(valid_data.size)

            if np.any(block_data.mask):
                hull_candidates = RasterDriver._hull_candidate_mask(~block_data.mask)
                hull_shapes = [geometry.shape(s) for s, _ in features.shapes(
                    np.ones(hull_candidates.shape, 'uint8'),
                    mask=hull_candidates,
                    transform=windows.transform(w, dataset.transform)
                )]
            else:
                w, s, e, n = windows.bounds(w, dataset.transform)
                hull_shapes = [geometry.Polygon([(w, s), (e, s), (e, n), (w, n)])]
            convex_hull = geometry.MultiPolygon([convex_hull, *hull_shapes]).convex_hull

            tdigest.update(valid_data)
            sstats.update(valid_data)

        if sstats.count() == 0:
            return None

        convex_hull_wgs = warp.transform_geom(
            dataset.crs, 'epsg:4326', geometry.mapping(convex_hull)
        )

        return {
            'valid_percentage': valid_data_count / total_count * 100,
            'range': (sstats.min(), sstats.max()),
            'mean': sstats.mean(),
            'stdev': sstats.std(),
            'percentiles': tdigest.quantile(np.arange(0.01, 1, 0.01)),
            'convex_hull': convex_hull_wgs
        } 

def _compute_image_stats(dataset: 'DatasetReader',
                             max_shape: Sequence[int] = None) -> Optional[Dict[str, Any]]:
        """Compute statistics for the given rasterio dataset by reading it into memory."""
        from rasterio import features, warp, transform
        from shapely import geometry

        out_shape = (dataset.height, dataset.width)

        if max_shape is not None:
            out_shape = (
                min(max_shape[0], out_shape[0]),
                min(max_shape[1], out_shape[1])
            )

        data_transform = transform.from_bounds(
            *dataset.bounds, height=out_shape[0], width=out_shape[1]
        )
        raster_data = dataset.read(1, out_shape=out_shape, masked=True)

        if dataset.nodata is not None:
            # nodata values might slip into output array if out_shape < dataset.shape
            raster_data[raster_data == dataset.nodata] = np.ma.masked

        valid_data = raster_data.compressed()

        if valid_data.size == 0:
            return None

        if np.any(raster_data.mask):
            hull_candidates = RasterDriver._hull_candidate_mask(~raster_data.mask)
            hull_shapes = (geometry.shape(s) for s, _ in features.shapes(
                np.ones(hull_candidates.shape, 'uint8'),
                mask=hull_candidates,
                transform=data_transform
            ))
            convex_hull = geometry.MultiPolygon(hull_shapes).convex_hull
        else:
            # no masked entries -> convex hull == dataset bounds
            w, s, e, n = dataset.bounds
            convex_hull = geometry.Polygon([(w, s), (e, s), (e, n), (w, n)])

        convex_hull_wgs = warp.transform_geom(
            dataset.crs, 'epsg:4326', geometry.mapping(convex_hull)
        )

        return {
            'valid_percentage': valid_data.size / raster_data.size * 100,
            'range': (float(valid_data.min()), float(valid_data.max())),
            'mean': float(valid_data.mean()),
            'stdev': float(valid_data.std()),
            'percentiles': np.percentile(valid_data, np.arange(1, 100)),
            'convex_hull': convex_hull_wgs
        }