Python rasterio.features() Examples

def __init__(self, raster_mask, no_data_value=None, ignore_labels=None):
        if ignore_labels is None:
            ignore_labels = []

        self.geometries = [{'label': int(label),
                            'polygon': Polygon(LinearRing(shp['coordinates'][0]),
                                               [LinearRing(pts) for pts in shp['coordinates'][1:]])}
                           for index, (shp, label) in enumerate(rasterio.features.shapes(raster_mask, mask=None))
                           if (int(label) is not no_data_value) and (int(label) not in ignore_labels)]

        self.areas = np.asarray([entry['polygon'].area for entry in self.geometries])
        self.decomposition = [shapely.ops.triangulate(entry['polygon']) for entry in self.geometries]

        self.label2cc = collections.defaultdict(list)
        for index, entry in enumerate(self.geometries):
            self.label2cc[entry['label']].append(index) 

def clip_rasters(folder_in, folder_out, aoi_in, debug=False):
    """Read continental rasters one at a time, clip to AOI and save

    Parameters
    ----------
    folder_in : str, Path
        Path to directory containing rasters.
    folder_out : str, Path
        Path to directory to save clipped rasters.
    aoi_in : str, Path
        Path to an AOI file (readable by Fiona) to use for clipping.
    """

    if isinstance(aoi_in, gpd.GeoDataFrame):
        aoi = aoi_in
    else:
        aoi = gpd.read_file(aoi_in)

    coords = [json.loads(aoi.to_json())["features"][0]["geometry"]]

    for file_path in os.listdir(folder_in):
        if file_path.endswith(".tif"):
            if debug:
                print(f"Doing {file_path}")
            ntl_rd = rasterio.open(os.path.join(folder_in, file_path))
            ntl, affine = mask(dataset=ntl_rd, shapes=coords, crop=True, nodata=0)

            if ntl.ndim == 3:
                ntl = ntl[0]

            save_raster(folder_out / file_path, ntl, affine) 

def rasterize_overlapped(self, shapes, out, **rasterize_args):
        """ Rasterize overlapped classes.

        :param shapes: Shapes to be rasterized.
        :type shapes: an iterable of pairs (rasterio.polygon, int)
        :param out: A numpy array to which to rasterize polygon classes.
        :type out: numpy.ndarray
        :param rasterize_args: Keyword arguments to be passed to `rasterio.features.rasterize`.
        :type rasterize_args: dict
        """
        rasters = [rasterio.features.rasterize([shape], out=np.copy(out), **rasterize_args) for shape in shapes]

        overlap_mask = np.zeros(out.shape, dtype=np.bool)
        no_data = self.no_data_value

        out[:] = rasters[0][:]
        for raster in rasters[1:]:
            overlap_mask[(out != no_data) & (raster != no_data) & (raster != out)] = True
            out[raster != no_data] = raster[raster != no_data]

        out[overlap_mask] = self.overlap_value 

def projectShapes(features, toCRS):
    import pyproj
    from functools import partial
    import fiona.crs as fcrs
    from shapely.geometry import shape, mapping
    from shapely.ops import transform as shpTrans

    project = partial(
        pyproj.transform,
        pyproj.Proj(fcrs.from_epsg(4326)),
        pyproj.Proj(toCRS))

    return list(
        {'geometry': mapping(
            shpTrans(
                project,
                shape(feat['geometry']))
        )} for feat in features
        ) 

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 

def test_compute_metadata_unoptimized(unoptimized_raster_file):
    from terracotta import exceptions
    from terracotta.drivers.raster_base import RasterDriver

    with rasterio.open(str(unoptimized_raster_file)) as src:
        data = src.read(1, masked=True)
        valid_data = data.compressed()
        dataset_shape = list(rasterio.features.dataset_features(
            src, bidx=1, band=False, as_mask=True, geographic=True
        ))

    convex_hull = MultiPolygon([shape(s['geometry']) for s in dataset_shape]).convex_hull

    # compare
    with pytest.warns(exceptions.PerformanceWarning):
        mtd = RasterDriver.compute_metadata(str(unoptimized_raster_file), use_chunks=False)

    np.testing.assert_allclose(mtd['valid_percentage'], 100 * valid_data.size / data.size)
    np.testing.assert_allclose(mtd['range'], (valid_data.min(), valid_data.max()))
    np.testing.assert_allclose(mtd['mean'], valid_data.mean())
    np.testing.assert_allclose(mtd['stdev'], valid_data.std())

    # allow some error margin since we only compute approximate quantiles
    np.testing.assert_allclose(
        mtd['percentiles'],
        np.percentile(valid_data, np.arange(1, 100)),
        rtol=2e-2
    )

    assert geometry_mismatch(shape(mtd['convex_hull']), convex_hull) < 1e-8 

def test_compute_metadata_nocrick(big_raster_file_nodata, monkeypatch):
    with rasterio.open(str(big_raster_file_nodata)) as src:
        data = src.read(1, masked=True)
        valid_data = data.compressed()
        dataset_shape = list(rasterio.features.dataset_features(
            src, bidx=1, band=False, as_mask=True, geographic=True
        ))

    convex_hull = MultiPolygon([shape(s['geometry']) for s in dataset_shape]).convex_hull

    from terracotta import exceptions
    import terracotta.drivers.raster_base

    with monkeypatch.context() as m:
        m.setattr(terracotta.drivers.raster_base, 'has_crick', False)

        with pytest.warns(exceptions.PerformanceWarning):
            mtd = terracotta.drivers.raster_base.RasterDriver.compute_metadata(
                str(big_raster_file_nodata), use_chunks=True
            )

    # compare
    np.testing.assert_allclose(mtd['valid_percentage'], 100 * valid_data.size / data.size)
    np.testing.assert_allclose(mtd['range'], (valid_data.min(), valid_data.max()))
    np.testing.assert_allclose(mtd['mean'], valid_data.mean())
    np.testing.assert_allclose(mtd['stdev'], valid_data.std())

    # allow error of 1%, since we only compute approximate quantiles
    np.testing.assert_allclose(
        mtd['percentiles'],
        np.percentile(valid_data, np.arange(1, 100)),
        rtol=2e-2
    )

    assert geometry_mismatch(shape(mtd['convex_hull']), convex_hull) < 1e-8 

def test_compute_metadata_approximate(nodata_type, big_raster_file_nodata, big_raster_file_mask):
    from terracotta.drivers.raster_base import RasterDriver

    if nodata_type == 'nodata':
        raster_file = big_raster_file_nodata
    elif nodata_type == 'masked':
        raster_file = big_raster_file_mask

    with rasterio.open(str(raster_file)) as src:
        data = src.read(1, masked=True)
        valid_data = data.compressed()
        dataset_shape = list(rasterio.features.dataset_features(
            src, bidx=1, band=False, as_mask=True, geographic=True
        ))

    convex_hull = MultiPolygon([shape(s['geometry']) for s in dataset_shape]).convex_hull

    # compare
    mtd = RasterDriver.compute_metadata(str(raster_file), max_shape=(512, 512))

    np.testing.assert_allclose(mtd['valid_percentage'], 100 * valid_data.size / data.size, atol=1)
    np.testing.assert_allclose(
        mtd['range'], (valid_data.min(), valid_data.max()), atol=valid_data.max() / 100
    )
    np.testing.assert_allclose(mtd['mean'], valid_data.mean(), rtol=0.02)
    np.testing.assert_allclose(mtd['stdev'], valid_data.std(), rtol=0.02)

    np.testing.assert_allclose(
        mtd['percentiles'],
        np.percentile(valid_data, np.arange(1, 100)),
        atol=valid_data.max() / 100, rtol=0.02
    )

    assert geometry_mismatch(shape(mtd['convex_hull']), convex_hull) < 0.05 

def __init__(self, n_samples, ref_mask_feature, ref_labels, sample_features, return_new_eopatch=False,
                 **sampling_params):
        """ Initialise sampling task.

        The data to be sampled is supposed to be a time-series stored in `DATA` type of the eopatch, while the raster
        image is supposed to be stored in `MASK_TIMELESS`. The output sampled features are stored in `DATA` and have
        shape T x N_SAMPLES x 1 x D, where T is the number of time-frames, N_SAMPLES the number of random samples, and D
        is the number of channels of the input time-series.

        The row and column index of sampled points can also be stored in the eopatch, to allow the same random sampling
        of other masks.

        :param n_samples: Number of random spatial points to be sampled from the time-series
        :type n_samples: int
        :param ref_mask_feature: Name of `MASK_TIMELESS` raster image to be used as a reference for sampling
        :type ref_mask_feature: str
        :param ref_labels: List of labels of `ref_mask_feature` mask which will be sampled
        :type ref_labels: list(int)
        :param sample_features: A collection of features that will be resampled. Each feature is represented by a tuple
            in a form of `(FeatureType, 'feature_name')` or
            (FeatureType, '<feature_name>', '<sampled feature name>'). If `sampled_feature_name` is not
            set the default name `'<feature_name>_SAMPLED'` will be used.

            Example: [(FeatureType.DATA, 'NDVI'), (FeatureType.MASK, 'cloud_mask', 'cloud_mask_1')]

        :type sample_features: list(tuple(FeatureType, str, str) or tuple(FeatureType, str))
        :param return_new_eopatch: If `True` the task will create new EOPatch, put sampled data and copy of timestamps
            and meta_info data in it and return it. If `False` it will just add sampled data to input EOPatch and
            return it.
        :type return_new_eopatch: bool
        :param sampling_params: Any other parameter used by `PointRasterSampler` class
        """
        self.n_samples = n_samples
        self.ref_mask_feature = self._parse_features(ref_mask_feature, default_feature_type=FeatureType.MASK_TIMELESS)
        self.ref_labels = list(ref_labels)
        self.sample_features = self._parse_features(sample_features, new_names=True,
                                                    rename_function='{}_SAMPLED'.format)
        self.return_new_eopatch = return_new_eopatch
        self.sampling_params = sampling_params 

def fillFacets(geoJSONpath, rasterPath, noProject, output, bands, zooming, batchprint, outputGeom, color):
    geoJSON, uidMap, featDims = getGJSONinfo(geoJSONpath)

    rasCRS, rasBounds, rasBands = getRasterInfo(rasterPath)

    bands = handleBandArgs(bands, rasBands)

    if rasCRS['proj'] == 'longlat' or noProject:
        noProject = True
        bounds = getBounds(geoJSON)
    else:
        ogeoJson = geoJSON
        geoJSON = projectShapes(geoJSON, rasCRS)
        bounds = getBounds(geoJSON)

    rasArr, oaff = loadRaster(rasterPath, bands, bounds)

    if min(rasArr.shape[0:2]) < 2 * featDims or zooming:
        rasArr = upsampleRaster(rasArr, featDims, zooming)


    if noProject:
        sampleVals = getRasterValues(geoJSON, rasArr, uidMap, bounds, outputGeom, bands, color)
    else:
        sampleVals = getRasterValues(geoJSON, rasArr, uidMap, bounds, outputGeom, bands, color, ogeoJson)

    if batchprint and outputGeom != True:
        sampleVals = batchStride(sampleVals, int(batchprint))

    if output:
        with open(output, 'w') as oFile:
            oFile.write(json.dumps({
                "type": "FeatureCollection",
                "features": list(sampleVals)
                }))
    else:
        try:
            for feat in sampleVals:
                click.echo(json.dumps(feat).rstrip())
        except IOError as e:
            pass 

def getGJSONinfo(geoJSONinfo):
    """
    Loads a lattice of GeoJSON, bounds, and creates a list mapping an on-the-fly UID w/ the actual index value.
    """
    features = list(i for i in filterBadJSON(geoJSONinfo))
    UIDs = list(feat['properties']['qt'] for feat in features)

    featDimensions = int(np.sqrt(len(features)/2.0))
    
    return features, UIDs, featDimensions 

def getBounds(features):
    xy = np.vstack(list(f['geometry']['coordinates'][0] for f in features))
    return coords.BoundingBox(
        xy[:,0].min(),
        xy[:,1].min(),
        xy[:,0].max(),
        xy[:,1].max()
        ) 

def execute(self, eopatch, seed=None):
        """ Execute random spatial sampling of time-series stored in the input eopatch

        :param eopatch: Input eopatch to be sampled
        :type eopatch: EOPatch
        :param seed: Setting seed of random sampling. If None no seed will be used.
        :type seed: int or None
        :return: An EOPatch with spatially sampled temporal features and associated labels
        :type eopatch: EOPatch
        """
        np.random.seed(seed)

        # Retrieve data and raster label image from eopatch
        f_type, f_name = next(self.ref_mask_feature(eopatch))
        raster = eopatch[f_type][f_name]

        # Initialise sampler
        sampler = PointRasterSampler(self.ref_labels, **self.sampling_params)
        # Perform sampling
        rows, cols = sampler.sample(raster, n_samples=self.n_samples)

        if self.return_new_eopatch:
            new_eopatch = EOPatch()
            new_eopatch.timestamp = eopatch.timestamp[:]
            new_eopatch.bbox = eopatch.bbox  # Should be copied
            new_eopatch.meta_info = eopatch.meta_info.copy()  # Should be deep copied - implement this in core
        else:
            new_eopatch = eopatch

        # Add sampled features
        for feature_type, feature_name, new_feature_name in self.sample_features(eopatch):

            if feature_type.is_time_dependent():
                sampled_data = eopatch[feature_type][feature_name][:, rows, cols, :]
            else:
                sampled_data = eopatch[feature_type][feature_name][rows, cols, :]
            new_eopatch[feature_type][new_feature_name] = sampled_data[..., np.newaxis, :]

        return new_eopatch 

def get_geometry_mask(width: int, height: int,
                      geometry: GeometryLike,
                      x_min: float, y_min: float, res: float) -> np.ndarray:
    geometry = convert_geometry(geometry)
    # noinspection PyTypeChecker
    transform = affine.Affine(res, 0.0, x_min,
                              0.0, -res, y_min + res * height)
    return rasterio.features.geometry_mask([geometry],
                                           out_shape=(height, width),
                                           transform=transform,
                                           all_touched=True,
                                           invert=True) 

def __init__(self, features, *, values=None, values_column='VALUE', raster_dtype=None, **rasterio_params):
        """
        :param features: One or more raster mask features which will be vectorized together with an optional new name
        of vector feature. If no new name is given the same name will be used.

        Examples:
            features=(FeatureType.MASK, 'CLOUD_MASK', 'VECTOR_CLOUD_MASK')

            features=[(FeatureType.MASK_TIMELESS, 'CLASSIFICATION'), (FeatureType.MASK, 'MONOTEMPORAL_CLASSIFICATION')]

        :type features: object supported by eolearn.core.utilities.FeatureParser class
        :param values: List of values which will be vectorized. By default is set to ``None`` and all values will be
            vectorized
        :type values: list(int) or None
        :param values_column: Name of the column in vector feature where raster values will be written
        :type values_column: str
        :param raster_dtype: If raster feature mask is of type which is not supported by ``rasterio.features.shapes``
            (e.g. ``numpy.int64``) this parameter is used to cast the mask into a different type
            (``numpy.int16``, ``numpy.int32``, ``numpy.uint8``, ``numpy.uint16`` or ``numpy.float32``). By default
            value of the parameter is ``None`` and no casting is done.
        :type raster_dtype: numpy.dtype or None
        :param: rasterio_params: Additional parameters to be passed to `rasterio.features.shapes`. Currently
            available is parameter `connectivity`.
        """
        self.feature_gen = self._parse_features(features, new_names=True)
        self.values = values
        self.values_column = values_column
        self.raster_dtype = raster_dtype
        self.rasterio_params = rasterio_params

        for feature_type, _, _ in self.feature_gen:
            if not (feature_type.is_spatial() and feature_type.is_discrete()):
                raise ValueError('Input features should be a spatial mask, but {} found'.format(feature_type)) 

def extract_features_and_labels(dataset, tile_size, only_cache=False):
    """For each satellite image and its corresponding shapefiles in the dataset create
    tiled features and labels."""
    features = []
    labels = []

    for geotiff_path, shapefile_paths in dataset:
        tiled_features, tiled_labels = create_tiled_features_and_labels(
            geotiff_path, shapefile_paths, tile_size, only_cache)

        features += tiled_features
        labels += tiled_labels

    return features, labels 

def preprocess_data(tile_size, dataset, only_cache=False):
    """Create features and labels for a given dataset. The features are tiles which contain
    the three RGB bands of the satellite image, so they have the form (tile_size, tile_size, 3).
    Labels are bitmaps with 1 indicating that the corresponding pixel in the satellite image
    represents water."""

    print('_' * 100)
    print("Start preprocessing data.")

    features_train, labels_train = extract_features_and_labels(
        dataset["train"], tile_size, only_cache)
    features_test, labels_test = extract_features_and_labels(
        dataset["test"], tile_size, only_cache)

    return features_train, features_test, labels_train, labels_test 

def merge_rasters(folder, percentile=70):
    """Merge a set of monthly rasters keeping the nth percentile value.

    Used to remove transient features from time-series data.

    Parameters
    ----------
    folder : str, Path
        Folder containing rasters to be merged.
    percentile : int, optional (default 70.)
        Percentile value to use when merging using np.nanpercentile.
        Lower values will result in lower values/brightness.

    Returns
    -------
    raster_merged : numpy array
        The merged array.
    affine : affine.Affine
        The affine transformation for the merged raster.
    """

    affine = None
    rasters = []

    for file in os.listdir(folder):
        if file.endswith(".tif"):
            ntl_rd = rasterio.open(os.path.join(folder, file))
            rasters.append(ntl_rd.read(1))

            if not affine:
                affine = ntl_rd.transform

    raster_arr = np.array(rasters)

    raster_merged = np.percentile(raster_arr, percentile, axis=0)

    return raster_merged, affine 

def polygonize(self, data=None, mask=None, connectivity=4, transform=None):
        """
        Yield (polygon, value) for each set of adjacent pixels of the same value.
        Wrapper around rasterio.features.shapes

        From rasterio documentation:

        Parameters
        ----------
        data : numpy ndarray
        mask : numpy ndarray
               Values of False or 0 will be excluded from feature generation.
        connectivity : 4 or 8 (int)
                       Use 4 or 8 pixel connectivity.
        transform : affine.Affine
                    Transformation from pixel coordinates of `image` to the
                    coordinate system of the input `shapes`.
        """
        if not _HAS_RASTERIO:
            raise ImportError('Requires rasterio module')
        if data is None:
            data = self.mask.astype(np.uint8)
        if mask is None:
            mask = self.mask
        if transform is None:
            transform = self.affine
        shapes = rasterio.features.shapes(data, mask=mask, connectivity=connectivity,
                                          transform=transform)
        return shapes 

def _vectorize_single_raster(self, raster, affine_transform, crs, timestamp=None):
        """ Vectorizes a data slice of a single time component

        :param raster: Numpy array or shape (height, width, channels)
        :type raster: numpy.ndarray
        :param affine_transform: Object holding a transform vector (i.e. geographical location vector) of the raster
        :type affine_transform: affine.Affine
        :param crs: Coordinate reference system
        :type crs: sentinelhub.CRS
        :param timestamp: Time of the data slice
        :type timestamp: datetime.datetime
        :return: Vectorized data
        :rtype: geopandas.GeoDataFrame
        """
        mask = None
        if self.values:
            mask = np.zeros(raster.shape, dtype=np.bool)
            for value in self.values:
                mask[raster == value] = True

        geo_list = []
        value_list = []
        for idx in range(raster.shape[-1]):
            for geojson, value in rasterio.features.shapes(raster[..., idx],
                                                           mask=None if mask is None else mask[..., idx],
                                                           transform=affine_transform, **self.rasterio_params):
                geo_list.append(shapely.geometry.shape(geojson))
                value_list.append(value)

        series_dict = {
            self.values_column: pd.Series(value_list),
            'geometry': GeoSeries(geo_list)
        }
        if timestamp is not None:
            series_dict['TIMESTAMP'] = pd.Series([timestamp] * len(geo_list))

        vector_data = GeoDataFrame(series_dict, crs=crs.pyproj_crs())

        if not vector_data.geometry.is_valid.all():
            vector_data.geometry = vector_data.geometry.buffer(0)

        return vector_data 

def __init__(self, vector_input, raster_feature, *, values=None, values_column=None, raster_shape=None,
                 raster_resolution=None, raster_dtype=np.uint8, no_data_value=0, write_to_existing=False,
                 overlap_value=None, buffer=0, **rasterio_params):
        """
        :param vector_input: Vector data to be used for rasterization. It can be given as a feature in `EOPatch` or
            as an independent geopandas `GeoDataFrame`.
        :type vector_input: (FeatureType, str) or GeoDataFrame
        :param raster_feature: New or existing raster feature into which data will be written. If existing raster
            raster feature is given it will by default take existing values and write over them.
        :type raster_feature: (FeatureType, str)
        :param values: If `values_column` parameter is specified then only polygons which have one of these specified
            values in `values_column` will be rasterized. It can be also left to `None`. If `values_column` parameter
            is not specified `values` parameter has to be a single number into which everything will be rasterized.
        :type values: list(int or float) or int or float or None
        :param values_column: A column in gived dataframe where values, into which polygons will be rasterized,
            are stored. If it is left to `None` then `values` parameter should be a single number into which
            everything will be rasterized.
        :type values_column: str or None
        :param raster_shape: Can be a tuple in form of (height, width) or an existing feature from which the spatial
            dimensions will be taken. Parameter `raster_resolution` can be specified instead of `raster_shape`.
        :type raster_shape: (int, int) or (FeatureType, str) or None
        :param raster_resolution: Resolution of raster into which data will be rasterized. Has to be given as a
            number or a tuple of numbers in meters. Parameter `raster_shape` can be specified instead or
            `raster_resolution`.
        :type raster_resolution: float or (float, float) or None
        :param raster_dtype: Data type of the obtained raster array, default is `numpy.uint8`.
        :type raster_dtype: numpy.dtype
        :param no_data_value: Default value of all other raster pixels into which no value will be rasterized.
        :type no_data_value: int or float
        :param write_to_existing: If `True` it will write to existing raster array and overwrite parts of its values.
            If `False` it will create a new raster array and remove the old one. Default is `False`.
        :param overlap_value: A value to override parts of raster where polygons of different classes overlap. If None,
            rasterization overlays polygon as it is the default behavior of `rasterio.features.rasterize`.
        :type overlap_value: raster's dtype
        :type write_to_existing: bool
        :param buffer: Buffer value passed to vector_data.buffer() before rasterization. If 0, no buffering is done.
        :type buffer: float
        :param: rasterio_params: Additional parameters to be passed to `rasterio.features.rasterize`. Currently
            available parameters are `all_touched` and `merge_alg`
        """

        self.vector_input, self.raster_feature = self._parse_main_params(vector_input, raster_feature)

        self.values = values
        self.values_column = values_column
        if values_column is None and (values is None or not isinstance(values, (int, float))):
            raise ValueError("One of parameters 'values' and 'values_column' is missing")

        self.raster_shape = raster_shape
        self.raster_resolution = raster_resolution
        if (raster_shape is None) == (raster_resolution is None):
            raise ValueError("Exactly one of parameters 'raster_shape' and 'raster_resolution' has to be specified")

        self.raster_dtype = raster_dtype
        self.no_data_value = no_data_value
        self.write_to_existing = write_to_existing
        self.rasterio_params = rasterio_params
        self.overlap_value = overlap_value
        self.buffer = buffer 

def test_compute_metadata(big_raster_file_nodata, big_raster_file_nomask,
                          nodata_type, use_chunks):
    from terracotta.drivers.raster_base import RasterDriver

    if nodata_type == 'nodata':
        raster_file = big_raster_file_nodata
    elif nodata_type == 'nomask':
        raster_file = big_raster_file_nomask

    if use_chunks:
        pytest.importorskip('crick')

    with rasterio.open(str(raster_file)) as src:
        data = src.read(1, masked=True)
        valid_data = data.compressed()
        dataset_shape = list(rasterio.features.dataset_features(
            src, bidx=1, band=False, as_mask=True, geographic=True
        ))

    convex_hull = MultiPolygon([shape(s['geometry']) for s in dataset_shape]).convex_hull

    # compare
    if nodata_type == 'nomask':
        with pytest.warns(UserWarning) as record:
            mtd = RasterDriver.compute_metadata(str(raster_file), use_chunks=use_chunks)
            assert 'does not have a valid nodata value' in str(record[0].message)
    else:
        mtd = RasterDriver.compute_metadata(str(raster_file), use_chunks=use_chunks)

    np.testing.assert_allclose(mtd['valid_percentage'], 100 * valid_data.size / data.size)
    np.testing.assert_allclose(mtd['range'], (valid_data.min(), valid_data.max()))
    np.testing.assert_allclose(mtd['mean'], valid_data.mean())
    np.testing.assert_allclose(mtd['stdev'], valid_data.std())

    # allow some error margin since we only compute approximate quantiles
    np.testing.assert_allclose(
        mtd['percentiles'],
        np.percentile(valid_data, np.arange(1, 100)),
        rtol=2e-2, atol=valid_data.max() / 100
    )

    assert geometry_mismatch(shape(mtd['convex_hull']), convex_hull) < 1e-8 

def mask_by_shape(infile, outfile, shapefile, to_db=False, datatype='float32',
                  rescale=True, min_value=0.000001, max_value=1, ndv=None,
                  description=True):

    # import shapefile geometries
    with fiona.open(shapefile, 'r') as file:
        features = [feature['geometry'] for feature in file
                    if feature['geometry']]

    # import raster
    with rasterio.open(infile) as src:
        out_image, out_transform = rasterio.mask.mask(src, features, crop=True)
        out_meta = src.meta.copy()
        out_image = np.ma.masked_where(out_image == ndv, out_image)

    
    # unmask array
    out_image = out_image.data
    
    # if to decibel should be applied
    if to_db is True:
        out_image = convert_to_db(out_image)

    if rescale:
        # if we scale to another d
        if datatype != 'float32':

            if datatype == 'uint8':
                out_image = scale_to_int(out_image, min_value, max_value, 'uint8')
            elif datatype == 'uint16':
                out_image = scale_to_int(out_image, min_value, max_value, 'uint16')

    out_meta.update({'driver': 'GTiff', 'height': out_image.shape[1],
                     'width': out_image.shape[2], 'transform': out_transform,
                     'nodata': ndv, 'dtype': datatype, 'tiled': True,
                     'blockxsize': 128, 'blockysize': 128})

    with rasterio.open(outfile, 'w', **out_meta) as dest:
        dest.write(out_image)
        if description:
            dest.update_tags(1, 
                    BAND_NAME='{}'.format(os.path.basename(infile)[:-4]))
            dest.set_band_description(1, 
                    '{}'.format(os.path.basename(infile)[:-4])) 

def sieve(
        self,
        size=2,
        mask=None,
        connectivity=4,
        file_path=None,
        driver="GTiff",
        nodata=None,
        dtype=None,
    ):
        """Replace pixels with their largest neighbor

        Thin wrapper around the rasterio.features.sieve method.

        Parameters
        ----------
        size : integer (default 2)
            Minimum number of contigous pixels to retain
        
        mask : ndarray (optional, default None)
            Values of False or 0 will be excluded from the sieving process
        
        connectivity : integer (default 4)
            Use 4 or 8 pixel connectivity for grouping pixels into features.
            Default is 4.

        file_path : str (optional, default None)
            Optional path to save calculated Raster object. If not
            specified then a tempfile is used.

        driver : str (default 'GTiff')
            Named of GDAL-supported driver for file export.

        nodata : any number (optional, default None)
            Nodata value for new dataset. If not specified then a nodata value is set
            based on the minimum permissible value of the Raster's data type.

        dtype : str (optional, default None)
            Optionally specify a numpy compatible data type when saving to file. If not
            specified, a data type is set based on the data type of the RasterLayer.

        Returns
        -------
        pyspatialml.RasterLayer
            Filled RasterLayer
        """
        arr = sieve(
            source=self.read(masked=True),
            size=size,
            mask=mask,
            connectivity=connectivity,
        )

        layer = self._write(arr, file_path, driver, dtype, nodata)

        return layer 

def create_bitmap(raster_dataset, shapefile_paths, satellite_path):
    """Create the bitmap for a given satellite image."""

    satellite_img_name = get_file_name(satellite_path)
    cache_file_name = "{}_water.tif".format(satellite_img_name)
    cache_path = os.path.join(WATER_BITMAPS_DIR, cache_file_name)
    try:
        # Try loading the water bitmap from cache.
        print("Load water bitmap from {}".format(cache_path))
        bitmap = load_bitmap(cache_path)
        bitmap[bitmap == 255] = 1
        return bitmap
    except IOError as e:
        print("No cache file found.")

    water_features = np.empty((0, ))

    print("Create bitmap for water features.")
    for shapefile_path in shapefile_paths:
        try:
            print("Load shapefile {}.".format(shapefile_path))
            with fiona.open(shapefile_path) as shapefile:
                # Each feature in the shapefile also contains meta information such as
                # wether the features is a lake or a river. We only care about the geometry
                # of the feature i.e. where it is located and what shape it has.
                geometries = [feature['geometry'] for feature in shapefile]

                water_features = np.concatenate(
                    (water_features, geometries), axis=0)
        except IOError as e:
            print("No shapefile found.")
            sys.exit(1)

    # Now that we have the vector data of all water features in our satellite image
    # we "burn it" into a new raster so that we get a B/W image with water features
    # in white and the rest in black. We choose the value 255 so that there is a stark
    # contrast between water and non-water pixels. This is only for visualisation
    # purposes. For the classifier we use 0s and 1s.
    bitmap_image = rasterio.features.rasterize(
        ((g, 255) for g in water_features),
        out_shape=raster_dataset.shape,
        transform=raster_dataset.transform)

    save_bitmap(cache_path, bitmap_image, raster_dataset)

    bitmap_image[bitmap_image == 255] = 1
    return bitmap_image 

def create_tiled_features_and_labels(geotiff_path,
                                     shapefile_paths,
                                     tile_size,
                                     only_cache=False):
    """Create the features and labels for a given satellite image and its shapefiles."""

    # Try to load tiles from cache.
    satellite_img_name = get_file_name(geotiff_path)
    cache_file_name = "{}_{}.pickle".format(satellite_img_name, tile_size)
    cache_path = os.path.join(TILES_DIR, cache_file_name)
    try:
        print("Load tiles from {}.".format(cache_path))
        with open(cache_path) as f:
            tiles = pickle.load(f)

        return tiles["features"], tiles["labels"]
    except IOError as e:
        if only_cache:
            raise
        print("Cache not available. Compute tiles.")

    # The provided satellite images have a different coordinate reference system as
    # the familiar WGS 84 which uses Latitude and Longitude. So we need to reproject
    # the satellite image to the WGS 84 coordinate reference system.
    dataset, wgs84_path = reproject_dataset(geotiff_path)
    bands = np.dstack(dataset.read())

    # For the given satellite image create a bitmap which has 1 at every pixel which corresponds
    # to water in the satellite image. In order to do this we use water polygons from OpenStreetMap.
    # The water polygons are stored in forms of shapefiles and are given by "shapefile_paths".
    water_bitmap = create_bitmap(dataset, shapefile_paths, geotiff_path)

    # Tile the RGB bands of the satellite image and the bitmap.
    tiled_bands = create_tiles(bands, tile_size, wgs84_path)
    tiled_bitmap = create_tiles(water_bitmap, tile_size, wgs84_path)

    # Due to the projection the satellite image in the GeoTIFF is not a perfect rectangle and the
    # remaining space on the edges is blacked out. When we overlay the GeoTIFF with the
    # shapefile it also overlays features for the blacked out parts which means that if we don't
    # remove these tiles the classifier will be fed with non-empty labels for empty features.
    tiled_bands, tiled_bitmap = remove_edge_tiles(tiled_bands, tiled_bitmap,
                                                  tile_size, dataset.shape)

    save_tiles(cache_path, tiled_bands, tiled_bitmap)

    return tiled_bands, tiled_bitmap 

def rasterize(self, shapes, out_shape=None, fill=0, out=None, transform=None,
                  all_touched=False, default_value=1, dtype=None):
        """
        Return an image array with input geometries burned in.
        Wrapper around rasterio.features.rasterize

        From rasterio documentation:

        Parameters
        ----------
        shapes : iterable of (geometry, value) pairs or iterable over
                 geometries.
        out_shape : tuple or list
                    Shape of output numpy ndarray.
        fill : int or float, optional
               Fill value for all areas not covered by input geometries.
        out : numpy ndarray
              Array of same shape and data type as `image` in which to store
              results.
        transform : affine.Affine
                    Transformation from pixel coordinates of `image` to the
                    coordinate system of the input `shapes`.
        all_touched : boolean, optional
                      If True, all pixels touched by geometries will be burned in.  If
                      false, only pixels whose center is within the polygon or that
                      are selected by Bresenham's line algorithm will be burned in.
        default_value : int or float, optional
                        Used as value for all geometries, if not provided in `shapes`.
        dtype : numpy data type
                Used as data type for results, if `out` is not provided.
        """
        if not _HAS_RASTERIO:
            raise ImportError('Requires rasterio module')
        if out_shape is None:
            out_shape = self.shape
        if transform is None:
            transform = self.affine
        raster = rasterio.features.rasterize(shapes, out_shape=out_shape, fill=fill,
                                             out=out, transform=transform,
                                             all_touched=all_touched,
                                             default_value=default_value, dtype=dtype)
        return raster