Mosaic with harmo 'band' : unexpected result

OTB Usage
1

Hi

My first post !

Context

I’m trying to make a mosaic with Mosaic Application and harmonization method to ‘band’.

Configuration setup

My system: ubuntu 20.04.1
Version of the OTB: 7.1
I installed the OTB with: don’t know…
Python version: 3.8.10
i’m using also pytob version: 1.4.0

Description of my issue

I’m doing a mosaic with three tiles of Sentinel-2. It’s not with raw bands but with a DHI (dynamic habitat index that is an annual synthesis of radiometric indices) of NDVI.

  • The mosaic with S2 raw bands and parameter 'harmo.method': 'band' is correct

  • The mosaic with DHIs and and paramter 'harmo.method': 'band', the mosaic is strange and stats are :
    MAX: 1056, MEAN: 223, MIN: -6844

Stats of the DHI used are:
tile T32TLR : MAX: 999, MEAN: 688, MIN: -493
tile T31TGL : MAX: 999, MEAN: 798, MIN: -998
tile T31TGK : MAX: 999, MEAN: 727, MIN: -994

Here is a capture of the three independent DHIs:

dhi_min_tiles
dhi_min_tiles790×953 153 KB

And the mosaic:

dhi_min_mosaic
dhi_min_mosaic819×935 52.5 KB

The gains seems not to be normal but I don’t understand why:
Gains : 0.625725 0.292308 0

For another combination of DHIs the gains are different:
Gains : 0 0 0.295846

  • The mosaic with DHIs and paramter 'harmo.method': 'none', the mosaic is correct but with no harmonisation… and at the connection between tiles there is no data values.

dhi_min_mosaic_none
dhi_min_mosaic_none820×933 152 KB

Here part of my script:

dhi_mosaic = pyotb.Mosaic({'il': tile_list,
                           'comp.feather': 'none',
                           'harmo.method': 'band',
                           'harmo.cost': 'rmse',
                           'interpolator' : 'bco',
                           'output.spacingx' : 100,
                           'output.spacingy' : 100,
                           'nodata': -10000
                           })
filename_extension = "?&streaming:type=tiled&streaming:sizemode=height&streaming:sizevalue=512"
dhi_mosaic.write(dhi_mosaic_path, pixel_type='int16', filename_extension=filename_extension)

Is someone can help me ?

Thanks !

2

Hi @Josselin

I don’t know DHI but I suspect that since it comes from NDVI, it should already be normalized. So not using any harmonization method at all is what I would recommend.

Anyway, I suspect that the wrong harmonization is because you decimate the image from its original resolution (10 meters?), creating some aliasing artifacts.
Since the rmse is used to compute the cost function, statistics are not very relevant since they are computed from images full of artifacts (using mu would have been a bit better since less sensible to pixel-wise spreads).

If you want something at a different scale, I encourage you to resample the images using some approach that estimates locally the value (e.g. average prior to decimation).
For instance say your DHI mosaic has 10m spacing, and you want some output at 100m spacing, you can do the following:

imgs = [pyotb.Smoothing({"in": i, "type": "mean", "radius": 5}) for i in tile_list]  # Average over 11x11 pixels everywhere (roughly 100m)
mosaic = pyotb.Mosaic({"il": tile_list, "nodata": -10000, "output.spacingx" : 100, "output.spacingy" : 100})

(After this you might want to post-process your output since Smoothing has polluted the boundaries of valid pixels with the -10000. Just process the binary mask separately, propagating boundaries using binary morpho box radius 5, mosaic, then re-apply it at the end).

If you really want to apply some harmonization I bet that resampling images prior to mosaicing them “works” (even if I am not convinced with the need of harmonization over normalized indice :angel: )

3

Thanks first

The DHI : CES Variables pour la biodiversité – Theia

Certainly the best answer to my problem because I need it at 10m spatial resolution.

Basically, I decimate from original spatial resolution in order to speed up the output writing and to enhance visualization on QGis for test…

So Smoothing is suitable than RigidTransformResample ?

4

Hi @Josselin ,

I never have used RigidTransformResample, so I can’t tell.
I don’t know if it could replace smoothing + decimation…