Blocks with value 0 in TSA output

[spengler-iamo]

When producing a TSS using a binary mask, blocks with value 0 are appearing in the output. Maybe related to #234

Parameterization:

++PARAM_TSA_START++

# INPUT/OUTPUT DIRECTORIES
# ------------------------------------------------------------------------
DIR_LOWER = /data/Aldhani/users/spengles/ReGreenUA/raster/force/level2
DIR_HIGHER = /data/Aldhani/users/spengles/ReGreenUA/raster/force/tss
DIR_PROVENANCE = /data/Aldhani/users/spengles/ReGreenUA/raster/force/provenance

# MASKING
# ------------------------------------------------------------------------
DIR_MASK = /data/Aldhani/users/spengles/ReGreenUA/raster/force/dc
BASE_MASK = shelterbelt_mask_3035.tif

# OUTPUT OPTIONS
# ------------------------------------------------------------------------
OUTPUT_FORMAT = GTiff
FILE_OUTPUT_OPTIONS = NULL
OUTPUT_EXPLODE = TRUE
OUTPUT_SUBFOLDERS = FALSE
FAIL_IF_EMPTY = FALSE

# PARALLEL PROCESSING
# ------------------------------------------------------------------------
NTHREAD_READ = 60
NTHREAD_COMPUTE = 20
NTHREAD_WRITE = 60
STREAMING = TRUE
PRETTY_PROGRESS = TRUE

# PROCESSING EXTENT AND RESOLUTION
# ------------------------------------------------------------------------
X_TILE_RANGE = 110 140
Y_TILE_RANGE = 10 80
FILE_TILE = NULL
CHUNK_SIZE = 1000 1000
RESOLUTION = 10
REDUCE_PSF = FALSE
USE_L2_IMPROPHE = FALSE

# SENSOR ALLOW-LIST
# ------------------------------------------------------------------------
SENSORS = SEN2A SEN2B SEN2C
TARGET_SENSOR = SEN2L
PRODUCT_TYPE_MAIN = BOA
PRODUCT_TYPE_QUALITY = QAI
SPECTRAL_ADJUST = FALSE

# QAI SCREENING
# ------------------------------------------------------------------------
SCREEN_QAI = NODATA CLOUD_OPAQUE CLOUD_BUFFER CLOUD_CIRRUS CLOUD_SHADOW SNOW SUBZERO SATURATION
ABOVE_NOISE = 0
BELOW_NOISE = 0

# PROCESSING TIMEFRAME
# ------------------------------------------------------------------------
DATE_RANGE = 2021-03-01 2025-10-31
DOY_RANGE = 1 365
DATE_IGNORE_LANDSAT_7 = 2099-12-31

# SPECTRAL INDEX
# ------------------------------------------------------------------------
INDEX = BLUE GREEN RED NIR SWIR1 SWIR2 REDEDGE1 REDEDGE2 REDEDGE3 BROADNIR
STANDARDIZE_TSS = NONE
OUTPUT_TSS = TRUE

# SPECTRAL MIXTURE ANALYSIS
# ------------------------------------------------------------------------
# This block only applies if INDEX includes SMA
# ------------------------------------------------------------------------
FILE_ENDMEM  = NULL
SMA_SUM_TO_ONE = TRUE
SMA_NON_NEG = TRUE
SMA_SHD_NORM = TRUE
SMA_ENDMEMBER = 1
OUTPUT_RMS = FALSE

# INTERPOLATION PARAMETERS
# ------------------------------------------------------------------------
INTERPOLATE = NONE
MOVING_MAX = 16
RBF_SIGMA = 8 16 32
RBF_CUTOFF = 0.95
HARMONIC_TREND = TRUE
HARMONIC_MODES = 3
HARMONIC_FIT_RANGE = 2015-01-01 2017-12-31
OUTPUT_NRT = FALSE
INT_DAY = 16
STANDARDIZE_TSI = NONE
OUTPUT_TSI = FALSE

# PYTHON UDF PARAMETERS
# ------------------------------------------------------------------------
FILE_PYTHON = NULL
PYTHON_TYPE = PIXEL
OUTPUT_PYP = FALSE

# R UDF PARAMETERS
# ------------------------------------------------------------------------
# Note: due to OpenMP threading conflicts, the usage of R UDFs will trigger
# the internal disabling of the streaming functionality (STREAMING = FALSE)
# ------------------------------------------------------------------------
FILE_RSTATS = NULL
RSTATS_TYPE = PIXEL
OUTPUT_RSP = FALSE

# SPECTRAL TEMPORAL METRICS
# ------------------------------------------------------------------------
OUTPUT_STM = FALSE
STM = Q25 Q50 Q75 AVG STD

# FOLDING PARAMETERS
# ------------------------------------------------------------------------
FOLD_TYPE = AVG
STANDARDIZE_FOLD = NONE
OUTPUT_FBY = FALSE
OUTPUT_FBQ = FALSE
OUTPUT_FBM = FALSE
OUTPUT_FBW = FALSE
OUTPUT_FBD = FALSE
OUTPUT_TRY = FALSE
OUTPUT_TRQ = FALSE
OUTPUT_TRM = FALSE
OUTPUT_TRW = FALSE
OUTPUT_TRD = FALSE
OUTPUT_CAY = FALSE
OUTPUT_CAQ = FALSE
OUTPUT_CAM = FALSE
OUTPUT_CAW = FALSE
OUTPUT_CAD = FALSE

# LAND SURFACE PHENOLOGY PARAMETERS - POLAR-BASED
# ------------------------------------------------------------------------
# Please note that these only work well with non-negative indices.--------
# ------------------------------------------------------------------------
POL_START_THRESHOLD = 0.2
POL_MID_THRESHOLD = 0.5
POL_END_THRESHOLD = 0.8
POL_ADAPTIVE = TRUE
POL = VSS VPS VES VSA RMR IGS
STANDARDIZE_POL = NONE
OUTPUT_PCT = FALSE
OUTPUT_POL = FALSE
OUTPUT_TRO = FALSE
OUTPUT_CAO = FALSE

# TREND PARAMETERS
# ------------------------------------------------------------------------
TREND_TAIL = TWO
TREND_CONF = 0.95
CHANGE_PENALTY = FALSE

++PARAM_TSA_END++

I am using FORCE v. 3.10.04 in Docker on Ubuntu 22.04.5 LTS

[davidfrantz]

Hi,

so far I cannot reproduce this.

• Can you please upload your mask for one tile where this occurs?
• Can you please post the FULL parameter file (also use double back-ticks to properly display the code block)
• Can you test the latest dev version?

[spengler-iamo]

I added the full parameter file (properly formatted, sorry for that)

Here is the mask for tile X0123_Y0047
shelterbelt_mask_3035.tif

Using v. 3.10.04-dev:::2026-04-02_11:02:47 I get the same output

[davidfrantz]

When you run this multiple times, are the block artefacts always in the same place?

[davidfrantz]

I got something... It might be some nasty effect in the combination of sparse Tiff files with the internal block layout.

[davidfrantz]

The shape of the black rectangle is one 256 x 256 px Block, i.e., 2560 m. In this region, I am seeing 0-values instead of proper NODATA values. The chunking size you chose for processing is smaller than that, i.e., 1 km.

[davidfrantz]

Can you confirm that your block effects are aligned with the 2560 x 2560 m block spacing?

[davidfrantz]

And can you please check whether this artifact occurs multiple times per image, or just once?

[spengler-iamo]

When I run it with same parameterization, the artefacts are always in the same place. Also all images in a tile contain the artefacts.

In this example I can confirm that the block effects are aligned with the 2560 x 2560 m spacing.

Initially I used 30000 x 30000 as chunk size. I also experienced the block artifacts with this chunk size, but at different locations compared to this example. I assumed that the artifacts were related to this warning: "/raster/force/level2/X0121_Y0037/20240709_LEVEL2_SEN2A_BOA.tif, band 1: Access window out of range in RasterIO(). " which appeared for multiple images.
I tried to overcome this by reducing the chunk size.

However I now tried with chunk size = 30000 x 30000 again and now I see no artefacts.

For troubleshooting I also repeated the cubing of the mask. Here I tried to change the internal block size of the mask before force-cube. The internal block size of the cubed mask files is now Block=3000x256.
Do I understand correctly that the chunk size must be higher than the internal block size of the mask?

I will redo the analysis now at broader scale and report back, if the artifacts are still appearing.

[spengler-iamo]

Now there are no artifacts in the output, when using CHUNK_SIZE = 30000 30000

[davidfrantz]

The underlying issue seems to be rooted in how GDAL deals with its internal blocks when writing sparse data.

Definitions (to be on the same page)

block: physical block structure within the Tiff (defaults to 256 x 256 px). Only the output Tiff matters in this case.

chunk: the piece of data that is used for processing.

Generally, the file layout is sparse. This means that blocks (the default 256px units), that are never touched by a processing chunk, are never written. This generally works fine, and the reader interprets those blocks as nodata.

The issue emerges when a chunk touches a block that was not written before. If the chunk-block-alignment is unfavorable, the piece of data covering the chunk is correctly updated in the Tiff, but the part of the block not covered by the chunk is set to 0, i.e., not honouring the nodata value.

I might need to refactor the block/chunk mechanic, unfortunately. It really is a pity that GTiff only supports multiples of 16 in their block sizes for tiled layouts. Strangely, this rule does not exist for striped layouts...

For now, setting the CHUNK_SIZE to TILE_SIZE is a viable workaround as it makes sure that the whole image (containing all blocks) is written.