Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sentinel2 L1C data wrongly scaled #1

Open
eomasters-repos opened this issue Aug 12, 2024 · 5 comments
Open

Sentinel2 L1C data wrongly scaled #1

eomasters-repos opened this issue Aug 12, 2024 · 5 comments

Comments

@eomasters-repos
Copy link

eomasters-repos commented Aug 12, 2024

Hello,

I think the Sentinel2-L1C data is wrongly scaled in the read_as_numpy function..

The np_arr contains the DN values and a scaled offset (default=-0.1) should not be applied.
I've noticed this when using sentinel2_superresolution and the values read where too low by 1000 (the scaling is set to 1).

Also, the no-data of the raw values is zero, but it is set to -10000 by default.

See page 444: https://sentinels.copernicus.eu/documents/247904/685211/S2-PDGS-TAS-DI-PSD-V14.9.pdf/3d3b6c9c-4334-dcc4-3aa7-f7c0deffbaf7?t=1643013091529

def read_as_numpy(
self,
bands: List[Band],
masks: Optional[List[Mask]] = None,
scale: float = 10000,
crs: Optional[str] = None,
resolution: float = 10,
no_data_value: float = np.nan,
bounds: Optional[rio.coords.BoundingBox] = None,
algorithm=rio.enums.Resampling.cubic,
dtype: np.dtype = np.dtype("float32"),
) -> Tuple[np.ndarray, Optional[np.ndarray], np.ndarray, np.ndarray, str,]:
"""Read bands from Sentinel2 products as a numpy
ndarray. Depending on the parameters, an internal WarpedVRT
dataset might be used.
:param bands: The list of bands to read
:param scale: Scale factor applied to reflectances (r_s = r /
scale). No scaling if set to None
:param crs: Projection in which to read the image (will use WarpedVRT)
:param resolution: Resolution of data. If different from the
resolution of selected bands, will use WarpedVRT
:param region: The region to read as a BoundingBox object or a
list of pixel coords (xmin, ymin, xmax, ymax)
:param no_data_value: How no-data will appear in output ndarray
:param bounds: New bounds for datasets. If different from
image bands, will use a WarpedVRT
:param algorithm: The resampling algorithm to be used if WarpedVRT
:param dtype: dtype of the output Tensor
:return: The image pixels as a np.ndarray of shape [bands,
width, height],
The masks pixels as a np.ndarray of shape [masks, width,
height],
The x coords as a np.ndarray of shape [width],
the y coords as a np.ndarray of shape [height],
the crs as a string
"""
if masks is None:
masks = self.ALL_MASKS
if len(bands):
img_files = [self.build_band_path(b) for b in bands]
np_arr, xcoords, ycoords, out_crs = regulargrid.read_as_numpy(
img_files,
crs=crs,
resolution=resolution,
offsets=self.offsets,
output_no_data_value=no_data_value,
input_no_data_value=-10000,
bounds=bounds,
algorithm=algorithm,
separate=True,
dtype=dtype,
scale=scale,
)
# Skip first dimension
np_arr = np_arr[0, ...]
# Apply radiometric offset
np_arr = np_arr + (self.radiometric_offset / scale)

@jmichel-otb
Copy link
Collaborator

I agree on the no_data which should be set to 0 on line 321. I do not agree on removing the offset however. The purpose of sensorsio is to provide ready to use surface reflectance by default, not DN. According to the documentation you pointed out (p. 444), RADIO_ADD_OFFSET should be added if available in metadata (logic handled on line 118).

If you would like to be able to disable this behaviour we can consider an extra flag in read_as_numpy(), and a corresponding switch in sentinel2_superresolution.

@eomasters-repos
Copy link
Author

I mean the offset should not be removed and also not the scaling.
But the scaling at Line#33 should not be:
np_arr = np_arr + (self.radiometric_offset / scale)
This is what is written in the PSD. But the PSD is wrong.
The formula results for a DN value of 5000 in 4999.9.
See here: https://sentiwiki.copernicus.eu/web/s2-products#S2Products-L1CS2-Products-L1Ctrue
Also discussed here: https://forum.step.esa.int/t/info-introduction-of-additional-radiometric-offset-in-pb04-00-products/35431/16?u=marco_eom

The line should be:
np_arr = (np_arr + self.radiometric_offset) / scale
which will give a value of 0.4

@jmichel-otb
Copy link
Collaborator

Ah, this is nasty. Thank you for spotting this. It will be fixed among other things in L1C driver (see https://framagit.org/jmichel-otb/sensorsio/-/merge_requests/14 ).

@jmichel-otb
Copy link
Collaborator

Fixes have been pushed to mr 14, branch is merged and github repo is synced.

Let me know if that works, so that I can close this issue.

@eomasters-repos
Copy link
Author

I haven't checkout the code. have only used the lib via sentinel2_superresolution. So can't test it right now. But changes look good to me.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants