EOdal Mapper: Grids do not align when reprojection is necessary from one UTM Zone into another

[lukasValentin]

Description:

When using the Mapper to fetch data from Sentinel-2 tiles in different UTM zones, the data needs to be reprojected from one UTM zone into another. This works fine. However, the resulting grids do not necessarily fully match because of the reprojection routine not "snapping" to a reference grid.

E.g., when data from UTM zones 32 and 33 is read, and the majority of scenes is in zone 32, the scenes in zone 33 projection are re-projected into zone 32N. However, the mapper currently does not use a reference grid, i.e., the grid of the 32 scenes, so that there might be small offsets in the coordinates.

This does not matter as long as the data is not stacked, e.g., for time series analysis.

How to recreate the problem

The code snippet below fetches S2 data from zones 31 and 32N (western Switerland).

from datetime import datetime
from pathlib import Path
from shapely.geometry import box
from typing import List

from eodal.config import get_settings
from eodal.core.scene import SceneCollection
from eodal.core.sensors.sentinel2 import Sentinel2
from eodal.mapper.feature import Feature
from eodal.mapper.filter import Filter
from eodal.mapper.mapper import Mapper, MapperConfigs

Settings = get_settings()
# set to False to use a local data archive
Settings.USE_STAC = True

# user-inputs
# -------------------------- Collection -------------------------------
collection: str = 'sentinel2-msi'

# ------------------------- Time Range ---------------------------------
time_start: datetime = datetime(2022, 6, 10)  		# year, month, day (incl.)
time_end: datetime = datetime(2022, 6, 30)   		# year, month, day (incl.)

# ---------------------- Spatial Feature  ------------------------------
bbox = box(*[6.5738, 46.4565, 7.2628, 47.2190])
feature = Feature(
    name='lake_neuchatel',
    geometry=bbox,
    epsg=4326,
    attributes={}
)

# ------------------------- Metadata Filters ---------------------------
metadata_filters: List[Filter] = [
    Filter('cloudy_pixel_percentage', '<', 25),
    Filter('processing_level', '==', 'Level-2A')
]

# query the scenes available (no I/O of scenes, this only fetches metadata)
mapper_configs = MapperConfigs(
    collection=collection,
    time_start=time_start,
    time_end=time_end,
    feature=feature,
    metadata_filters=metadata_filters)

# to enhance reproducibility and provide proper documentation, the MapperConfigs
# can be saved as yaml (and also then be loaded again from yaml)
mapper_configs.to_yaml('data/sample_mapper_call.yaml')

# now, a new Mapper instance is created
mapper = Mapper(mapper_configs)
mapper.query_scenes()

# load the least cloudy scene available from STAC
scene_kwargs = {
    'scene_constructor': Sentinel2.from_safe,
    'scene_constructor_kwargs': {'band_selection': ["B02"], 'read_scl': False}}

mapper.load_scenes(scene_kwargs=scene_kwargs)

ulx_list = []
uly_list = []
epsg_list = []
for _, scene in mapper.data:
    ulx = scene['blue'].geo_info.ulx
    uly = scene['blue'].geo_info.uly
    epsg = scene['blue'].geo_info.epsg
    ulx_list.append(ulx)
    uly_list.append(uly)
    epsg_list.append(epsg)

# resulting ulx_list -> there is an offset of ~2 m to 9.6 m
>>> [770564.8761871913, 770564.8761871913, 770566.8953721962, 770564.8761871913, 770574.5245121085, 770574.5245121085]
# resulting uly_list -> larger offsets here are a due no-data in single tiles
>>> [5238323.835322518, 5238323.835322518, 5238322.759971117, 5238323.835322518, 5207434.29725682, 5207434.29725682]

As mentioned also by @orianif it would be nicer if all scenes in a SceneCollection share the same reference grid to allow stacks through time.

Proposed solution

We should think about a post-processing step in the Mapper to ensure all scenes are projected into the same reference grid. This will then also ensure that all scenes share exactly the same extent in terms of rows and columns.

[lukasValentin]

Update

I'll add a post-processing step to the Mapper class using Band.reproject to reproject all bands and hence scenes into a common reference grid and extent (i.e., common for all entries in the SceneCollection).

The common reference grid is determined by collecting the bounds of all scenes in the collection. The total_bound of all these bounds defines then the reference. We can then use Band.reproject using a destination raster with the extent of the common reference grid (determined from the maximum number of rows and columns within the common reference bounds). Interpolation if any (in case of reprojection from different CRS, e.g., different UTM zones) is done using nearest neighbor interpolation.

IMPORTANT NOTICE: It turned out that this will only work (as the mosaicing by the way) if all bands in the scenes have the same spatial resolution. Thus, you should either make sure to only query bands that have the same spatial resolution or use the scene_modifier custom-function to re-sample all bands into a common spatial resolution.

[lukasValentin]

closed by #65