added capabilities to calculate indices

[lucadelu]

added capabilities to calculate indices from level 2 output.

Right now the script is working but the result is all 0, so something is not working properly. I tried to calculate the same index using QGIS and it works properly.

@griembauer the gdal_calc.py call is the following one, do you see any problem on it?

gdal_calc.py -A /mnt/cri_digiagri/open/test/20230202_LEVEL2_SEN2A_BOA_clearsky.tif --a_band=4 -B 20230202_LEVEL2_SEN2A_BOA_clearsky.tif --b_band=8 --outfile=20230202_LEVEL2_SEN2A_NDVI.tif --calc="((A.astype(float)-B)/(A.astype(float)+B))" --type=Float32 --creation-option=COMPRESS=ZSTD --creation-option=TILED=YES --creation-option=PREDICTOR=2 --creation-option=BIGTIFF=YES

otherwise I was thinking to use rasterio to calculate the index

[griembauer]

Note that FORCE omits the ultra-blue, cirrus, and water vapor bands, so to compute the NDVI you would need to use bands 3 (red) and 7 (broad NIR). See the Sentinel-2 output table here: https://force-eo.readthedocs.io/en/latest/components/lower-level/level2/format.html

Apart from that, I see no issue with the command - I can also reproduce the issue that you are facing, with the output being 0 only, but I don't have an idea how to fix it yet.

[griembauer]

I think I found the issue: the --A_band and --B_band parameters must be with capital letters in the beginning (in your example it is --a_band and --b_band). Also it should be (B-A)/(B+A) if B is NIR and A is Red, but this has nothing to do with the issues being 0 only

This example worked for me:

gdal_calc.py -A 20210627_LEVEL2_SEN2B_BOA_clearsky.tif --A_band=3 -B 20210627_LEVEL2_SEN2B_BOA_clearsky.tif --B_band=7 --outfile=20210627_LEVEL2_SEN2B_NDVI.tif --calc="((B.astype(float)-A)/(B.astype(float)+A))" --type=Float32 --creation-option=COMPRESS=ZSTD --creation-option=TILED=YES --creation-option=PREDICTOR=2 --creation-option=BIGTIFF=YES