How to use the rasterio.features function in rasterio

def test_shapes():
    """Test creation of shapes from pixel values"""

    image = numpy.zeros((20, 20), dtype=rasterio.ubyte)
    image[5:15, 5:15] = 127
    with rasterio.drivers():
        shapes = ftrz.shapes(image)
        shape, val = next(shapes)
        assert shape['type'] == 'Polygon'
        assert len(shape['coordinates']) == 2  # exterior and hole
        assert val == 0
        shape, val = next(shapes)
        assert shape['type'] == 'Polygon'
        assert len(shape['coordinates']) == 1  # no hole
        assert val == 127
        try:
            shape, val = next(shapes)
        except StopIteration:
            assert True
        else:
            assert False

image = numpy.zeros((rows, cols), dtype=dtype)
                image[2:5, 2:5] = test_value
                ftrz.sieve(image, 2)

        # Test mask types
        image = numpy.zeros((rows, cols), dtype='uint8')
        image.fill(255)

        supported_mask_types = (
            ('bool', 1),
            ('uint8', 255)
        )
        for dtype, mask_value in supported_mask_types:
            mask = numpy.zeros((rows, cols), dtype=dtype)
            mask[2:5, 2:5] = mask_value
            sieved_image = ftrz.sieve(image, 2, mask=mask)
            assert numpy.array_equal(image, sieved_image)

        unsupported_mask_types = (
            ('int8', -127),
            ('int16', -32768)
        )
        for dtype, mask_value in unsupported_mask_types:
            with pytest.raises(ValueError):
                mask = numpy.zeros((rows, cols), dtype=dtype)
                mask[2:5, 2:5] = mask_value
                ftrz.sieve(image, 2, mask=mask)

if do_transform:
            affine_obj = reference_im.transform

    # extract geometries and pair them with burn values
    if burn_field:
        if out_type == 'int':
            feature_list = list(zip(buffered_df[geom_col],
                                    buffered_df[burn_field].astype('uint8')))
        else:
            feature_list = list(zip(buffered_df[geom_col],
                                    buffered_df[burn_field].astype('uint8')))
    else:
        feature_list = list(zip(buffered_df[geom_col],
                                [burn_value] * len(buffered_df)))

    output_arr = features.rasterize(shapes=feature_list, out_shape=shape,
                                    transform=affine_obj)
    if min_background_value:
        output_arr = np.clip(output_arr, min_background_value,
                             np.max(output_arr))

    if out_file:
        meta = reference_im.meta.copy()
        meta.update(count=1)
        if out_type == 'int':
            meta.update(dtype='uint8')
        with rasterio.open(out_file, 'w', **meta) as dst:
            dst.write(output_arr, indexes=1)

    return output_arr

for fname in images:
        # ensure formats match
        with rasterio.open(str(fname), 'r') as ds:
            transform = ds.transform
            img = ds.read(1)

            if mask_value is not None:
                new_mask = img & mask_value == mask_value
            else:
                new_mask = img != ds.nodata
            if mask is None:
                mask = new_mask
            else:
                mask |= new_mask

    shapes = rasterio.features.shapes(mask.astype('uint8'), mask=mask)
    shape = shapely.ops.unary_union([shapely.geometry.shape(shape) for shape, val in shapes if val == 1])

    # convex hull
    geom = shape.convex_hull

    # buffer by 1 pixel
    geom = geom.buffer(1, join_style=3, cap_style=3)

    # simplify with 1 pixel radius
    geom = geom.simplify(1)

    # intersect with image bounding box
    geom = geom.intersection(shapely.geometry.box(0, 0, mask.shape[1], mask.shape[0]))

    # transform from pixel space into CRS space
    geom = shapely.affinity.affine_transform(geom, (transform.a, transform.b, transform.d,

def geometry_mask(geoms, geobox, all_touched=False, invert=False):
    """
    Create a mask from shapes.

    By default, mask is intended for use as a
    numpy mask, where pixels that overlap shapes are False.
    :param list[Geometry] geoms: geometries to be rasterized
    :param datacube.utils.GeoBox geobox:
    :param bool all_touched: 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.
    :param bool invert: If True, mask will be True for pixels that overlap shapes.
    """
    return rasterio.features.geometry_mask([geom.to_crs(geobox.crs) for geom in geoms],
                                           out_shape=geobox.shape,
                                           transform=geobox.affine,
                                           all_touched=all_touched,
                                           invert=invert)

lst_vector = [vector for vector in src]

    # Sort feature in order to priorize the burning in the raster image (ex: vegetation before roads...)
    if attribute_name is not None:
        lst_vector.sort(key=lambda vector: get_key_recursive(attribute_name, vector))

    lst_vector_tuple = lst_ids(list_vector=lst_vector, attr_name=attribute_name, target_ids=target_ids, merge_all=merge_all)

    if merge_all:
        return rasterio.features.rasterize([v for vecs in lst_vector_tuple.values() for v in vecs],
                                           fill=fill,
                                           out_shape=input_image.shape,
                                           transform=input_image.transform,
                                           dtype=np.int16)
    else:
        burned_rasters = [rasterio.features.rasterize(lst_vector_tuple[id],
                                                      fill=fill,
                                                      out_shape=input_image.shape,
                                                      transform=input_image.transform,
                                                      dtype=np.int16) for id in lst_vector_tuple]
        return np.stack(burned_rasters, axis=-1)

if cartoCSS:
        for i in breaks:
            click.echo('[value = ' + str(breaks[i]) + '] { polygon-fill: @class' + str(i) + '}')

    schema = { 'geometry': 'MultiPolygon', 'properties': { 'value': 'float' } }

    with fiona.open(outfile, "w", "ESRI Shapefile", schema, crs=src.crs) as outshp:
        tRas = np.zeros(classRas.shape, dtype=np.uint8)
        click.echo("Vectorizing: ", nl=False)
        for i, br in enumerate(breaks):
            click.echo("%d, " % (br), nl=False)
            tRas[np.where(classRas>=i)] = 1
            tRas[np.where(classRas 5 or c == 0:
                            if axonometrize:
                                f = np.array(f)
                                f[:,1] += (axonometrize * br)
                            if nosimple:
                                 poly = Polygon(f)
                            else:
                                poly = Polygon(f).simplify(simplest / float(smoothing), preserve_topology=True)
                            featurelist.append(poly)
                    if len(featurelist) != 0:
                        oPoly = MultiPolygon(featurelist)
                        outshp.write({'geometry': mapping(oPoly),'properties': {'value': br}})

featureList = ((geom, int(value)) for geom, value in zip(srcGDF.geometry, srcGDF[burnValueField]))
    else:
        featureList = []

    with rasterio.open(srcRasterFileName) as rst:
        meta = rst.meta.copy()
        meta.update(count=1)
        meta.update(dtype='uint8')

        with rasterio.open(
                outRasterFileName, 'w',
                **meta) as dst:

            if featureList:
                print(featureList)
                burned = features.rasterize(shapes=featureList,
                                            out_shape=rst.shape,
                                            transform=rst.transform
                                           )
            else:
                burned = np.empty(rst.shape, dtype='uint8')

            dst.write(burned, indexes=1)

    return srcGDF

This only works for 1D lat and lon arrays.

        for internal use: does not check valitity of input
    """
    # TODO: use only this function once https://github.com/mapbox/rasterio/issues/1844
    # is resolved

    from rasterio import features

    shapes = zip(polygons, numbers)

    transform = _transform_from_latlon(lon, lat)
    out_shape = (len(lat), len(lon))

    raster = features.rasterize(
        shapes,
        out_shape=out_shape,
        fill=fill,
        transform=transform,
        dtype=np.float,
        **kwargs
    )

    return raster

tgdf.to_file(temp_shp_fn)        
        
        # Get the template raster
        rst = rasterio.open(rst_fn)

        # copy and update the metadata from the input raster for the output
        meta = rst.meta.copy()
        meta.update(compress='lzw')

        with rasterio.open(out_fn, 'w+', **meta) as out:
            out_arr = out.read(1)

            # this is where we create a generator of geom, value pairs to use in rasterizing
            shapes = ((geom,value) for geom, value in zip(tgdf.geometry, tgdf.RZ_key))

            burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=out.transform)
            out.write_band(1, burned)