OSGeo / grass-addons

GRASS GIS Addons Repository
https://grass.osgeo.org/grass-stable/manuals/addons/
GNU General Public License v2.0
104 stars 154 forks source link

[Feat] i.sentinel.import: handle processing baselines #1019

Open ninsbl opened 9 months ago

ninsbl commented 9 months ago

Name of the addon Processing baselines of Sentinel-2 are regularly updated. Some changes may pose significant changes in the data (as e.g. baseline 4.00): https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/processing-baseline

i.sentinel.import should issue a warning when data with new processing baselines are encountered that the module is not known to handle. So, a list of "supported Baselines" should be added and checked against. That increases the need for maintenance a bit, but avoids surprises for users...

As for the BOA offsett introduced in baseline 4.00, I would suggest that i.sentinel.import accounts for that (if requested by the user) and rescales accordingly...

Expected behavior Users are made aware of changes in processing baselines, and the module is updated/adjusted if needed.

ninsbl commented 9 months ago

Not sure if this actually should be flagged as a bug, as i.sentinel.preproc and i.sentinel.mask will likely fail with un-adjusted baseline 4.00 data...

ecodiv commented 5 months ago

@ninsbl Does i.sentinel.import do anything with the rescale (QUANTIFICATION_VALUE) factor? If I import Sentinel-2 data, the imported bands have the unscaled (original) value. And I can't find where rescaling should take place in the i.sentinel.import code.

The page you link to gives the following information for processing version 5.0: " The Processing Baseline identifier 05.00 tags every Sentinel-2 Collection-1 product issued from the reprocessing activity. As part of the Sentinel-2 reprocessing campaign, the improvements introduced in recent processing baselines 04.00 deployed on 25 January 2022 is generalised to the historical archive". That suggests that older images now also are based on the new processing baselines? Indeed, if I download an image from e.g., 2021, the MTD-MSIL2A.xml file does provide the BOA_QUANTIFICATION_VALUE unit=10000 and BOA_ADD_OFFSET band_id=-1000 (for all bands).

Do you know how this affects the use of i.sentinel.mask? For the latter, would it just be a matter of adding the option to provide an offset value, besides the existing option to provide the rescale factor? Or should the rules used in i.sentinel.mask be adjusted to account for possible negative values?

ninsbl commented 5 months ago

Yes, this affects i.sentinel.mask and all indices like NDWI, NDSI, NDWI, ...

This is how I rescale in a custom python script:

def rescale_boa(map_name, offset=-1000.0, quantization=10000.0, environment=None):
    """Rescale input to BOA
    From the user’s point of view, the L1C Top of Atmosphere (TOA) reflectance (L1C_TOA) shall be retrieved from the output radiometry as follows:
      - Digital Number DN=0 will remain the “NO_DATA” value
      - For a given DN in [1;215-1], the L1C TOA reflectance value will be:
        "L1C_TOAi = (L1C_DNi + RADIO_ADD_OFFSETi) / QUANTIFICATION_VALUEi"
    """
    expression = f"{TMP_NAME}_{map_name.split('@')[0]}=float(if({map_name}+{offset} < 0, 0.0, if({map_name}+{offset} > {quantization}, 1.0, float({map_name}+{offset}) / {quantization})))"
    gs.run_command(
        "r.mapcalc",
        expression=expression,
        overwrite=True,
        #quiet=True,
        env=environment,
    )
    return f"{TMP_NAME}_{map_name.split('@')[0]}"

Where _TMPNAME is a global variable...

Rescaling should probably be added here: https://github.com/OSGeo/grass-addons/blob/5ae27ab8023af3010475e3719cd516f12fd2281d/src/imagery/i.sentinel/i.sentinel.import/i.sentinel.import.py#L289-L320

That said, nowadays, I would probably use GDALs Sentinel-2 driver: https://gdal.org/drivers/raster/sentinel2.html for import and simplify i.sentinel.import.

Another option to consider here is to use GDAL VRTs to handle rescaling and possibly projection "on the fly" (meaning, before importing)...

In the Sentinel-2 driver description, there is a section on VRTs and you can read more here: https://gdal.org/drivers/raster/vrt.html

Here an example for creating (also warped (=reprojected)) VRT files: https://github.com/OSGeo/grass-addons/blob/5ae27ab8023af3010475e3719cd516f12fd2281d/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py#L657

The nice thing with VRT and the sentinel driver is that zip-files do not have to be extracted, You can either link or import and reproject in one go, It should be possible to even rescale, reproject and import in one step... And the Sentinel-2 driver reads the metadata too...