improve cutline handling

[yellowcap]

This is another attempt to resolve #588

The proposed changes ensure that the FirstMethod finishes when all pixels over the target geometry is filled.

The main implementation challenge is: the pixel_selection class needs to know the target mask to use in is_done function.

This means that the input geometry needs to be converted to a mask as numpy array to which the pixel_selection needs to have access.

For this case, the cutline vrt_option approach does not work, because for gdal cutlines the mask generation happens on the gdal side and the mask is not returned. So the proposed changes do not use the vrt_option argument, but rasterized the geometry directly and applies this mask on the python side after requesting pixels.

import datetime

import httpx
import matplotlib.pyplot as plt
from rasterio.crs import CRS
from rasterio.warp import transform_geom

from rio_tiler.io import STACReader
from rio_tiler.mosaic import mosaic_reader

stac_endpoint = "https://earth-search.aws.element84.com/v0/search"

# Date filter
date_min = "2022-08-01"
date_max = "2022-08-31"

start = datetime.datetime.strptime(date_min, "%Y-%m-%d").strftime("%Y-%m-%dT00:00:00Z")
end = datetime.datetime.strptime(date_max, "%Y-%m-%d").strftime("%Y-%m-%dT23:59:59Z")

# Bern large
aoi = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [296.138689305058506, -3217.936052156386722],
                [2770.964285295431182, -5479.41461400965818],
                [3496.34420136157496, -8466.273091929069778],
                [2045.584369229287404, -11794.486824467847327],
                [1106.857419026044226, -15506.725218453404523],
                [4435.071151564818138, -16701.46860962116989],
                [10110.10225961170363, -16232.105134519548301],
                [10408.78810740364861, -15336.047591143724276],
                [10238.110480093964725, -12903.891401980770752],
                [9811.416411819765926, -11367.792756193644891],
                [7763.284884103592049, -7271.529700761304412],
                [5629.814542732583504, -4113.993595532210747],
                [3795.030049153516302, -1980.523254161202203],
                [1490.882080472823873, -401.755201546653552],
                [296.138689305058506, -3217.936052156386722],
            ]
        ],
    },
}
epsg = 21780

query = {
    "collections": [
        "sentinel-s2-l2a-cogs"
    ],  # Make sure to query only sentinel-2 COGs collection
    "datetime": f"{start}/{end}",
    "query": {"eo:cloud_cover": {"lt": 90}},
    "intersects": transform_geom(f"epsg:{epsg}", "epsg:4326", aoi["geometry"]),
    "limit": 2,
    "fields": {
        "include": [
            "id",
            "properties.datetime",
            "properties.eo:cloud_cover",
        ],
        "exclude": ["links"],
    },
}

headers = {
    "Content-Type": "application/json",
    "Accept-Encoding": "gzip",
    "Accept": "application/geo+json",
}


data = httpx.post(stac_endpoint, headers=headers, json=query).json()
print(data["context"])

sceneid = [f["id"] for f in data["features"]]
cloudcover = [f["properties"]["eo:cloud_cover"] for f in data["features"]]
dates = [f["properties"]["datetime"][0:10] for f in data["features"]]

items = [feat | {"stac_version": data["stac_version"]} for feat in data["features"]]


def tiler(asset, *args, **kwargs):
    print(asset["id"])
    with STACReader(None, asset) as src:
        img = src.feature(
            shape=kwargs["shape"],
            shape_crs=kwargs["shape_crs"],
            assets=kwargs["assets"],
        )
        print("min/max/shape", img.data.min(), img.data.max(), img.data.shape)

        return img

img, assets_used = mosaic_reader(
    items,
    tiler,
    shape=aoi,
    shape_crs=CRS.from_epsg(21780),
    assets=["B04", "B03", "B02"],
    threads=1,
)

fig = plt.figure(figsize=(30, 10))

ax = fig.add_subplot(1, 2, 1)
img.rescale(in_range=((0, 5000),))
ax.imshow(img.data_as_image())

ax = fig.add_subplot(1, 2, 2)
ax.imshow(img.mask)

plt.show()

[yellowcap]

For the docs, it seems like aside from the all_touched part, there is no api change on the functions. So maybe just adding another example would be a good way to document?

Either an extension of https://cogeotiff.github.io/rio-tiler/examples/Using-rio-tiler-mosaic/

or a fresh one similar to the snipped in the ticket above.

[vincentsarago]

👍 all good for me.

Would you be able to write a note in the changelog 🙏

[yellowcap]

Added changelog note

[yellowcap]

One sidenote: the create_cutline utils function is no longer used in this codebase. What to do in this case? Just leave it be, or deprecate?

[vincentsarago]

@yellowcap why is the create_cutline no more used? I think it should still be!

[yellowcap]

@yellowcap why is the create_cutline no more used? I think it should still be!

I first left it in, but the problem is that the cutline mask is done in the projection of the source image, so it won't exactly align with the cutline mask done in the target projection. The "finish early" will fail in some of these cases because there are pixels that are masked in the original projection, but not anymore in the target one. So these pixels will always be masked in the target projection.

I thought about this in terms of what it means for data transfer efficiency. My guess is that this will not have implications for data transfers, but not sure. Maybe some COG internal tile requests can be skipped with the cutline vrt_option in place. Do you know if GDAL optimizes that under the hood, i.e. identifying what tiles not to get based on the cutline?

[yellowcap]

Also, I just realized that the masking will not fully work for rasters with multiple resolutions (10m vs 20m) in this case. So I'll have to digg into that, the stop early fails when multiple resolutions are included.

[yellowcap]

I added another commit that solves the multiple resolution problem. The approach is to only apply the same cutline mask to all assets using the highest resolution rasterization.
0894fad

[vincentsarago]

I thought about this in terms of what it means for data transfer efficiency. My guess is that this will not have implications for data transfers, but not sure. Maybe some COG internal tile requests can be skipped with the cutline vrt_option in place. Do you know if GDAL optimizes that under the hood, i.e. identifying what tiles not to get based on the cutline?

It might (really doubt) but I don't think this will be a major issueQ 👍

Let's keep create_cutline in utils because it might be useful outside the feature method (e.g to mask tiles developmentseed/titiler#376)

[vincentsarago]

🤔 what happens if there is no cutline_mask in the ImageData object? I think this line shouldn't be changed

[yellowcap]

You are absolutely right again! But that opens the problem that pixels might never get unmasked when upsampling. In the following example, the green pixel would always be masked at the coarse resolution, but at the higher resolution it would not be masked. So the finish early would not work again. The solution is to go back to what you already suggested earlier - setting always all_touched to True. In that case the finer mask should always be included in the coarser mask. Will make change accordingly.

[vincentsarago]

👍 we could also add a mechanism to set an option for finish early when like 95% of the cutline is filled 🤷

[vincentsarago]

👍

[vincentsarago]

@yellowcap all seems good to me now! do you want to add a tests for the multiple resolution case?

[yellowcap]

ok test for multires added, took the blue.tif test file, resampled from 900m to 1800m, and renamed the resampled file to lowres.tif. Added this as an asset to the stac.json item file, and using the lowres file in the test.

[vincentsarago]

We need to handle projection difference between the input feature and output crs

[vincentsarago]

thanks a ton @yellowcap for working on this 🙏

[yellowcap]

Awesome! Thanks for all the help and getting it over the last humps @vincentsarago !