Python rasterio.transform() Examples

def _reproject(self, eopatch, src_raster):
        """
        Reprojects the raster data from Geopedia's CRS (POP_WEB) to EOPatch's CRS.
        """
        height, width = src_raster.shape

        dst_raster = np.ones((height, width), dtype=self.raster_dtype)

        src_bbox = eopatch.bbox.transform(CRS.POP_WEB)
        src_transform = rasterio.transform.from_bounds(*src_bbox, width=width, height=height)

        dst_bbox = eopatch.bbox
        dst_transform = rasterio.transform.from_bounds(*dst_bbox, width=width, height=height)

        rasterio.warp.reproject(src_raster, dst_raster,
                                src_transform=src_transform, src_crs={'init': CRS.ogc_string(CRS.POP_WEB)},
                                src_nodata=0,
                                dst_transform=dst_transform, dst_crs={'init': CRS.ogc_string(eopatch.bbox.crs)},
                                dst_nodata=self.no_data_val)

        return dst_raster 

def bounds_to_ranges(out_bounds=None, in_affine=None, in_shape=None):
    """
    Return bounds range values from geolocated input.

    Parameters
    ----------
    out_bounds : tuple
        left, bottom, right, top
    in_affine : Affine
        input geolocation
    in_shape : tuple
        input shape

    Returns
    -------
    minrow, maxrow, mincol, maxcol
    """
    return itertools.chain(
        *from_bounds(
            *out_bounds, transform=in_affine, height=in_shape[-2], width=in_shape[-1]
        ).round_lengths(pixel_precision=0).round_offsets(pixel_precision=0).toranges()
    ) 

def __init__(
        self, in_tile=None, in_data=None, out_profile=None, out_tile=None, tags=None
    ):
        """Prepare data & profile."""
        out_tile = out_tile or in_tile
        validate_write_window_params(in_tile, out_tile, in_data, out_profile)
        self.data = extract_from_array(
            in_raster=in_data,
            in_affine=in_tile.affine,
            out_tile=out_tile
        )
        # use transform instead of affine
        if "affine" in out_profile:
            out_profile["transform"] = out_profile.pop("affine")
        self.profile = out_profile
        self.tags = tags 

def thin(guess_in):
    """
    Use scikit-image skeletonize to 'thin' the guess raster.

    Parameters
    ----------
    guess_in : path-like or 2D array
        Output from threshold().

    Returns
    -------
    guess_skel : numpy array
        Thinned version.
    affine : Affine
        Only if path-like supplied.
    """

    if isinstance(guess_in, (str, Path)):
        guess_rd = rasterio.open(guess_in)
        guess_arr = guess_rd.read(1)
        affine = guess_rd.transform

        guess_skel = skeletonize(guess_arr)
        guess_skel = guess_skel.astype("int32")

        return guess_skel, affine

    elif isinstance(guess_in, np.ndarray):
        guess_skel = skeletonize(guess_in)
        guess_skel = guess_skel.astype("int32")

        return guess_skel

    else:
        raise ValueError 

def _get_wms_request(self, bbox, size_x, size_y):
        """
        Returns WMS request.
        """
        bbox_3857 = bbox.transform(CRS.POP_WEB)

        return GeopediaWmsRequest(layer=self.layer,
                                  theme=self.theme,
                                  bbox=bbox_3857,
                                  width=size_x,
                                  height=size_y,
                                  image_format=self.image_format,
                                  custom_url_params={CustomUrlParam.TRANSPARENT: True}) 

def execute(self, eopatch):
        """ Execute function which adds new vector layer to the EOPatch

        :param eopatch: input EOPatch
        :type eopatch: EOPatch
        :return: New EOPatch with added vector layer
        :rtype: EOPatch
        """
        for raster_ft, raster_fn, vector_fn in self.feature_gen(eopatch):
            vector_ft = FeatureType.VECTOR_TIMELESS if raster_ft.is_timeless() else FeatureType.VECTOR

            raster = eopatch[raster_ft][raster_fn]
            height, width = raster.shape[:2] if raster_ft.is_timeless() else raster.shape[1: 3]

            if self.raster_dtype:
                raster = raster.astype(self.raster_dtype)

            affine_transform = rasterio.transform.from_bounds(*eopatch.bbox, width=width, height=height)

            crs = eopatch.bbox.crs

            if raster_ft.is_timeless():
                eopatch[vector_ft][vector_fn] = self._vectorize_single_raster(raster, affine_transform, crs)
            else:
                gpd_list = [self._vectorize_single_raster(raster[time_idx, ...], affine_transform, crs,
                                                          timestamp=eopatch.timestamp[time_idx])
                            for time_idx in range(raster.shape[0])]

                eopatch[vector_ft][vector_fn] = GeoDataFrame(pd.concat(gpd_list, ignore_index=True),
                                                             crs=gpd_list[0].crs)

        return eopatch 

def execute(self, eopatch):
        """ Execute method

        :param eopatch: input EOPatch
        :type eopatch: EOPatch
        :return: New EOPatch with vector data transformed into a raster feature
        :rtype: EOPatch
        """
        if eopatch.bbox is None:
            raise ValueError('EOPatch has to have a bounding box')

        height, width = self._get_raster_shape(eopatch)

        group_classes = self.overlap_value is not None
        rasterization_shapes = self._get_rasterization_shapes(eopatch, group_classes=group_classes)

        if not rasterization_shapes:
            eopatch[self.raster_feature] = np.full((height, width, 1), self.no_data_value, dtype=self.raster_dtype)
            return eopatch

        affine_transform = rasterio.transform.from_bounds(*eopatch.bbox, width=width, height=height)
        rasterize_args = dict(self.rasterio_params, transform=affine_transform, dtype=self.raster_dtype)

        raster = self._get_raster(eopatch, height, width)
        rasterize_func = rasterio.features.rasterize if self.overlap_value is None else self.rasterize_overlapped
        rasterize_func(rasterization_shapes, out=raster, **rasterize_args)

        eopatch[self.raster_feature] = raster[..., np.newaxis]
        return eopatch 

def pansharpen(vis, vis_transform, pan, pan_transform,
               pan_dtype, r_crs, dst_crs, weight,
               method="Brovey", src_nodata=0):
    """Pansharpen a lower-resolution visual band

    Parameters
    =========
    vis: ndarray, 3D with shape == (3, vh, vw)
        Visual band array with RGB bands
    vis_transform: Affine
        affine transform defining the georeferencing of the vis array
    pan: ndarray, 2D with shape == (ph, pw)
        Panchromatic band array
    pan_transform: Affine
        affine transform defining the georeferencing of the pan array
    method: string
        Algorithm for pansharpening; default Brovey

    Returns:
    ======
    pansharp: ndarray, 3D with shape == (3, ph, pw)
        pansharpened visual band
        affine transform is identical to `pan_transform`
    """
    rgb = _upsample(_create_apply_mask(vis), pan.shape, vis_transform, r_crs,
                    pan_transform, dst_crs)

    # Main Pansharpening Processing
    if method == "Brovey":
        pansharp, _ = Brovey(rgb, pan, weight, pan_dtype)
    # TODO: add other methods

    return pansharp 

def __init__(self, array, grid):
        from rasterio.transform import from_origin

        self.grid = grid
        self.cell_size = self.grid.cell_size

        self.array = array
        self._affine = from_origin(
            self.grid.cloud.data.min[0],
            self.grid.cloud.data.max[1],
            self.grid.cell_size,
            self.grid.cell_size,
        ) 

def __init__(self, path, cloud):
        import rasterio

        self.in_raster = rasterio.open(path)

        # Check cell size
        cell_size_x, cell_size_y = (
            self.in_raster.transform[0],
            abs(self.in_raster.transform[4]),
        )
        if cell_size_x != cell_size_y:
            print("Cell sizes not equal of input raster, not supported.")
            raise ValueError
        else:
            cell_size = cell_size_x

        self.cloud = cloud
        self.cell_size = cell_size

        min_x, max_x = self.in_raster.bounds[0], self.in_raster.bounds[2]
        min_y, max_y = self.in_raster.bounds[1], self.in_raster.bounds[3]

        self.m = self.in_raster.height
        self.n = self.in_raster.width

        # Create bins
        bins_x = np.searchsorted(
            np.linspace(min_x, max_x, self.n), self.cloud.data.points["x"]
        )
        bins_y = np.searchsorted(
            np.linspace(min_y, max_y, self.m), self.cloud.data.points["y"]
        )

        self.cloud.data.points["bins_x"] = bins_x
        self.cloud.data.points["bins_y"] = bins_y
        self.cells = self.cloud.data.points.groupby(["bins_x", "bins_y"]) 

def save_merged_rasters(self, datatype, fileroot=None):
        import rasterio
        from rasterio.merge import merge
        from rasterio.plot import show
        from os.path import join
        import argparse
        import glob
        
        if fileroot == None:
            fileroot = datatype

        root = '{}/{}/{}*.tif'.format(self.save_samples_path, datatype, fileroot)
        filename = '{}/{}/merged_{}.tif'.format(self.save_samples_path, datatype, fileroot)

        files = glob.glob(join(root))
        mosaic_rasters = [rasterio.open(file) for file in files]

        mosaic, out_transform = merge(mosaic_rasters)

        meta = (rasterio.open(files[0])).meta

        meta.update({"driver": "GTiff",
                    "height": mosaic.shape[1],
                    "width":  mosaic.shape[2],
                    "transform": out_transform})

        with rasterio.open(filename, "w", **meta) as dest:
            dest.write(mosaic)

        filename = '{}/{}/{}_merged.png'.format(self.save_samples_path, datatype, datatype) 
        self.save_raster_png(mosaic, filename)

        if 'output' in filename or 'target' in filename:
            self.save_height_colormap(filename, mosaic) 

def save_merged_rasters(self, datatype, fileroot=None):
        import rasterio
        from rasterio.merge import merge
        from rasterio.plot import show
        from os.path import join
        import argparse
        import glob
        
        if fileroot == None:
            fileroot = datatype

        root = '{}/{}/{}*.tif'.format(self.save_samples_path, datatype, fileroot)
        filename = '{}/{}/merged_{}.tif'.format(self.save_samples_path, datatype, fileroot)

        files = glob.glob(join(root))
        mosaic_rasters = [rasterio.open(file) for file in files]

        mosaic, out_transform = merge(mosaic_rasters)

        meta = (rasterio.open(files[0])).meta

        meta.update({"driver": "GTif",
                    "height": mosaic.shape[1],
                    "width":  mosaic.shape[2],
                    "transform": out_transform})

        with rasterio.open(filename, "w", **meta) as dest:
            dest.write(mosaic)

        filename = '{}/{}/{}_merged.png'.format(self.save_samples_path, datatype, datatype) 
        self.save_raster_png(mosaic, filename) 

def _check_alignment(layers):
        """Check that a list of raster datasets are aligned with the same pixel
        dimensions and geotransforms.

        Parameters
        ----------
        layers : list
            List of pyspatialml.RasterLayer objects.

        Returns
        -------
        dict or False
            Dict of metadata if all layers are spatially aligned, otherwise
            returns False.
        """

        src_meta = []
        for layer in layers:
            src_meta.append(layer.ds.meta.copy())

        if not all(i["crs"] == src_meta[0]["crs"] for i in src_meta):
            Warning(
                "crs of all rasters does not match, " "possible unintended consequences"
            )

        if not all(
            [
                i["height"] == src_meta[0]["height"]
                or i["width"] == src_meta[0]["width"]
                or i["transform"] == src_meta[0]["transform"]
                for i in src_meta
            ]
        ):
            return False

        else:
            return src_meta[0] 

def run(self, processes=4):
        """TODO"""
        if processes == 1:
            self.pool = MockTub(init_worker, (self.inpaths, self.global_args))
        else:
            self.pool = Pool(processes, init_worker, (self.inpaths, self.global_args))

        self.options["transform"] = guard_transform(self.options["transform"])

        if self.mode == "manual_read":
            reader_worker = manual_reader(self.run_function)
        elif self.mode == "array_read":
            reader_worker = array_reader(self.run_function)
        else:
            reader_worker = simple_reader(self.run_function)

        if isinstance(self.outpath_or_dataset, rasterio.io.DatasetWriter):
            destination = self.outpath_or_dataset
        else:
            destination = rasterio.open(self.outpath_or_dataset, "w", **self.options)

        # Open an output file, work through the function in parallel,
        # and write out the data.
        with destination as dst:
            for data, window in self.pool.imap_unordered(reader_worker, self.windows):
                dst.write(data, window=window)

        self.pool.close()
        self.pool.join() 

def threshold(dists_in, cutoff=0.0):
    """Convert distance array into binary array of connected locations.

    Parameters
    ----------
    dists_in : path-like or numpy array
        2D array output from gridfinder algorithm.
    cutoff : float, optional (default 0.5.)
        Cutoff value below which consider the cells to be grid.

    Returns
    -------
    guess : numpy array
        Binary representation of input array.
    affine: affine.Affine
        Affine transformation for raster.
    """
    if isinstance(dists_in, (str, Path)):
        dists_rd = rasterio.open(dists_in)
        dists_r = dists_rd.read(1)
        affine = dists_rd.transform

        guess = dists_r.copy()
        guess[dists_r > cutoff] = 0
        guess[dists_r <= cutoff] = 1

        return guess, affine

    elif isinstance(dists_in, np.ndarray):
        guess = dists_in.copy()
        guess[dists_in > cutoff] = 0
        guess[dists_in <= cutoff] = 1

        return guess

    else:
        raise ValueError 

def process_tile(tile):
    """Process a single MBTiles tile

    Parameters
    ----------
    tile : mercantile.Tile

    Returns
    -------

    tile : mercantile.Tile
        The input tile.
    bytes : bytearray
        Image bytes corresponding to the tile.

    """
    global base_kwds, resampling, src

    # Get the bounds of the tile.
    ulx, uly = mercantile.xy(
        *mercantile.ul(tile.x, tile.y, tile.z))
    lrx, lry = mercantile.xy(
        *mercantile.ul(tile.x + 1, tile.y + 1, tile.z))

    kwds = base_kwds.copy()
    kwds['transform'] = transform_from_bounds(ulx, lry, lrx, uly,
                                              kwds['width'], kwds['height'])
    src_nodata = kwds.pop('src_nodata', None)
    dst_nodata = kwds.pop('dst_nodata', None)

    warnings.simplefilter('ignore')

    with MemoryFile() as memfile:

        with memfile.open(**kwds) as tmp:

            # determine window of source raster corresponding to the tile
            # image, with small buffer at edges
            try:
                west, south, east, north = transform_bounds(TILES_CRS, src.crs, ulx, lry, lrx, uly)
                tile_window = window_from_bounds(west, south, east, north, transform=src.transform)
                adjusted_tile_window = Window(
                    tile_window.col_off - 1, tile_window.row_off - 1,
                    tile_window.width + 2, tile_window.height + 2)
                tile_window = adjusted_tile_window.round_offsets().round_shape()

                # if no data in window, skip processing the tile
                if not src.read_masks(1, window=tile_window).any():
                    return tile, None

            except ValueError:
                log.info("Tile %r will not be skipped, even if empty. This is harmless.", tile)

            reproject(rasterio.band(src, tmp.indexes),
                      rasterio.band(tmp, tmp.indexes),
                      src_nodata=src_nodata,
                      dst_nodata=dst_nodata,
                      num_threads=1,
                      resampling=resampling)

        return tile, memfile.read() 

def _rasterio_read(
    input_file=None,
    indexes=None,
    dst_bounds=None,
    dst_shape=None,
    dst_crs=None,
    resampling=None,
    src_nodata=None,
    dst_nodata=None,
):

    def _read(
        src, indexes, dst_bounds, dst_shape, dst_crs, resampling, src_nodata, dst_nodata
    ):
        height, width = dst_shape[-2:]
        if indexes is None:
            dst_shape = (len(src.indexes), height, width)
            indexes = list(src.indexes)
        src_nodata = src.nodata if src_nodata is None else src_nodata
        dst_nodata = src.nodata if dst_nodata is None else dst_nodata
        dst_left, dst_bottom, dst_right, dst_top = dst_bounds
        with WarpedVRT(
            src,
            crs=dst_crs,
            src_nodata=src_nodata,
            nodata=dst_nodata,
            width=width,
            height=height,
            transform=affine_from_bounds(
                dst_left, dst_bottom, dst_right, dst_top, width, height
            ),
            resampling=Resampling[resampling]
        ) as vrt:
            return vrt.read(
                window=vrt.window(*dst_bounds),
                out_shape=dst_shape,
                indexes=indexes,
                masked=True
            )
    if isinstance(input_file, str):
        logger.debug("got file path %s", input_file)
        try:
            with rasterio.open(input_file, "r") as src:
                return _read(
                    src, indexes, dst_bounds, dst_shape, dst_crs, resampling, src_nodata,
                    dst_nodata
                )
        except RasterioIOError as e:
            try:
                if path_exists(input_file):
                    raise e
            except:
                raise e
            raise FileNotFoundError("%s not found" % input_file)
    else:
        logger.debug("assuming file object %s", input_file)
        return _read(
            input_file, indexes, dst_bounds, dst_shape, dst_crs, resampling,
            src_nodata, dst_nodata
        ) 

def atmos(
    ctx,
    atmo,
    contrast,
    bias,
    jobs,
    out_dtype,
    src_path,
    dst_path,
    creation_options,
    as_color,
):
    """Atmospheric correction
    """
    if as_color:
        click.echo(
            "rio color {} {} {}".format(
                src_path, dst_path, simple_atmo_opstring(atmo, contrast, bias)
            )
        )
        exit(0)

    with rasterio.open(src_path) as src:
        opts = src.profile.copy()
        windows = [(window, ij) for ij, window in src.block_windows()]

    opts.update(**creation_options)
    opts["transform"] = guard_transform(opts["transform"])

    out_dtype = out_dtype if out_dtype else opts["dtype"]
    opts["dtype"] = out_dtype

    args = {"atmo": atmo, "contrast": contrast, "bias": bias, "out_dtype": out_dtype}

    jobs = check_jobs(jobs)

    if jobs > 1:
        with riomucho.RioMucho(
            [src_path],
            dst_path,
            atmos_worker,
            windows=windows,
            options=opts,
            global_args=args,
            mode="manual_read",
        ) as mucho:
            mucho.run(jobs)
    else:
        with rasterio.open(dst_path, "w", **opts) as dest:
            with rasterio.open(src_path) as src:
                rasters = [src]
                for window, ij in windows:
                    arr = atmos_worker(rasters, window, ij, args)
                    dest.write(arr, window=window) 

def rasterio_to_gdir(gdir, input_file, output_file_name,
                     resampling='cubic'):
    """Reprojects a file that rasterio can read into the glacier directory.

    Parameters
    ----------
    gdir : :py:class:`oggm.GlacierDirectory`
        the glacier directory
    input_file : str
        path to the file to reproject
    output_file_name : str
        name of the output file (must be in cfg.BASENAMES)
    resampling : str
        nearest', 'bilinear', 'cubic', 'cubic_spline', or one of
        https://rasterio.readthedocs.io/en/latest/topics/resampling.html
    """

    output_file = gdir.get_filepath(output_file_name)
    assert '.tif' in output_file, 'output_file should end with .tif'

    if not gdir.has_file('dem'):
        raise InvalidWorkflowError('Need a dem.tif file to reproject to')

    with rasterio.open(input_file) as src:

        kwargs = src.meta.copy()
        data = src.read(1)

        with rasterio.open(gdir.get_filepath('dem')) as tpl:

            kwargs.update({
                'crs': tpl.crs,
                'transform': tpl.transform,
                'width': tpl.width,
                'height': tpl.height
            })

            with rasterio.open(output_file, 'w', **kwargs) as dst:
                for i in range(1, src.count + 1):

                    dest = np.zeros(shape=(tpl.height, tpl.width),
                                    dtype=data.dtype)

                    reproject(
                        source=rasterio.band(src, i),
                        destination=dest,
                        src_transform=src.transform,
                        src_crs=src.crs,
                        dst_transform=tpl.transform,
                        dst_crs=tpl.crs,
                        resampling=getattr(Resampling, resampling)
                    )

                    dst.write(dest, indexes=i) 

def _polygon_to_pix(polygon):
    """Transforms polygon coordinates to integer pixel coordinates. It makes
    the geometry easier to handle and reduces the number of points.

    Parameters
    ----------
    polygon: the shapely.geometry.Polygon instance to transform.

    Returns
    -------
    a shapely.geometry.Polygon class instance.
    """

    def project(x, y):
        return np.rint(x).astype(np.int64), np.rint(y).astype(np.int64)

    poly_pix = shapely.ops.transform(project, polygon)

    # simple trick to correct invalid polys:
    tmp = poly_pix.buffer(0)

    # sometimes the glacier gets cut out in parts
    if tmp.type == 'MultiPolygon':
        # If only small arms are cut out, remove them
        area = np.array([_tmp.area for _tmp in tmp])
        _tokeep = np.argmax(area).item()
        tmp = tmp[_tokeep]

        # check that the other parts really are small,
        # otherwise replace tmp with something better
        area = area / area[_tokeep]
        for _a in area:
            if _a != 1 and _a > 0.05:
                # these are extremely thin glaciers
                # eg. RGI40-11.01381 RGI40-11.01697 params.d1 = 5. and d2 = 8.
                # make them bigger until its ok
                for b in np.arange(0., 1., 0.01):
                    tmp = shapely.ops.transform(project, polygon.buffer(b))
                    tmp = tmp.buffer(0)
                    if tmp.type == 'MultiPolygon':
                        continue
                    if tmp.is_valid:
                        break
                if b == 0.99:
                    raise InvalidGeometryError('This glacier geometry is not '
                                               'valid.')

    if not tmp.is_valid:
        raise InvalidGeometryError('This glacier geometry is not valid.')

    return tmp 

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 accuracy(grid_in, guess_in, aoi_in, buffer_amount=0.01):
    """Measure accuracy against a specified grid 'truth' file.

    Parameters
    ----------
    grid_in : str, Path
        Path to vector truth file.
    guess_in : str, Path
        Path to guess output from guess2geom.
    aoi_in : str, Path
        Path to AOI feature.
    buffer_amount : float, optional (default 0.01.)
        Leeway in decimal degrees in calculating equivalence.
        0.01 DD equals approximately 1 mile at the equator.
    """

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

    grid = gpd.read_file(grid_in)
    grid_clipped = clip_line_poly(grid, aoi)
    grid_buff = grid_clipped.buffer(buffer_amount)

    guesses_reader = rasterio.open(guess_in)
    guesses = guesses_reader.read(1)

    grid_for_raster = [(row.geometry) for _, row in grid_clipped.iterrows()]
    grid_raster = rasterize(
        grid_for_raster,
        out_shape=guesses_reader.shape,
        fill=1,
        default_value=0,
        all_touched=True,
        transform=guesses_reader.transform,
    )
    grid_buff_raster = rasterize(
        grid_buff,
        out_shape=guesses_reader.shape,
        fill=1,
        default_value=0,
        all_touched=True,
        transform=guesses_reader.transform,
    )

    grid_raster = flip_arr_values(grid_raster)
    grid_buff_raster = flip_arr_values(grid_buff_raster)

    tp = true_positives(guesses, grid_buff_raster)
    fn = false_negatives(guesses, grid_raster)

    return tp, fn