0

I need to rescale Sentinel2 data (both lvl 1C and 2A) to 8 bit, because of the requirements of a GRASS GIS module (i.histo.match) that only works with values from 0:255. All the documentation of Sentinel2 1C and 2A data products say that the data are encoded at 12bit, but that is not up to date (also this other post pointed this out https://gis.stackexchange.com/a/234044/172030). In fact, if I run gdalinfo on my L1C files, it says that it has 15 significant bits:

Image Structure Metadata:
COMPRESSION=JPEG2000
NBITS=15

while if I run it on my L2A files (created with Sen2Cor, not downloaded) it omits the "NBITS" field.

Image Structure Metadata:
COMPRESSION=JPEG2000

I don't know how to correctly perform the conversion: should I treat the data as formatted in 15 or 16 bits? I.e. DN / 2^15 * 2^8 or DN / 2^16 * 2^8? Is there any better way to perform this transformation?

Improve this question
2

You should be able to use the max= parameter to i.histo.match to retain the full range of values with input rasters like S2 that use 15 bit DN, without any rescaling.

I checked two S2 images (band 8 for two dates), using r.univar and found the max values about 16,000. So I ran:

micha@RMS:Western Negev$ i.histo.match input=`g.list rast sep=comma pat=b08*` max=18000 --o

The resulting histogram-matched rasters showed more or less the same range of values:

micha@RMS:Western Negev$ r.univar b08_20210125.match
 100%
total null and non-null cells: 49343168
total null cells: 44910908

Of the non-null cells:
----------------------
n: 4432260
minimum: 11
maximum: 16248
range: 16237
mean: 2356.4
mean of absolute values: 2356.4
standard deviation: 1625.54
variance: 2.64238e+06
variation coefficient: 68.9838 %
sum: 10444197995
Improve this answer
1
  • Thanks a lot. I don't know how I didn't see that, as I am usually obsessively careful to every small detail, but I guess this sometimes brings one not to see the forest for the trees! Thank you again – F.H. Jan 27 at 11:23
1

Here a script I have written some time ago which rescales to 8 bit and exports into GeoTIFF (or COG):

#!/bin/sh

# Markus Neteler, mundialis 2018
# 8bit export of satellite data, e.g. Sentinel-2

#BAND_R=T32ULB_20170314T104011_B04
#BAND_G=T32ULB_20170314T104011_B03
#BAND_B=T32ULB_20170314T104011_B02
BAND_R=mosaik_B08
BAND_G=mosaik_B03
BAND_B=mosaik_B02

# vector
AOI=ireland_aoi1
AOIEXT=gpkg

v.import input=$AOI.$ext out=$AOI

##### 
# rescale to 8bit
for band in $BAND_R $BAND_G $BAND_B ; do
  g.region vector=$AOI align=$band -p
  #eval `r.info $band -r`
  # leave 0 for nodata
  r.rescale input=$band output=${band}_8bit to=1,255
done

i.colors.enhance r=${BAND_R}_8bit g=${BAND_G}_8bit b=${BAND_B}_8bit

# we try to export with color table
i.group group=rgb_8bit input=${BAND_R}_8bit,${BAND_G}_8bit,${BAND_B}_8bit

# separate 8bit band export (optionally, write out COG)
for band in $BAND_R $BAND_G $BAND_B ; do
   r.out.gdal input=${band}_8bit output=${band}_geotiff_8bit.tif type=Byte create="COMPRESS=LZW" nodata=0 -m
done

# combine to RGB GeoTIFF
gdal_merge.py  -separate ${BAND_R}_geotiff_8bit.tif ${BAND_G}_geotiff_8bit.tif ${BAND_B}_geotiff_8bit.tif -o s2_${AOI}_geotiff_8bit.tif -co COMPRESS=LZW

# GeoTIFF -> TIFF (optionally, write out COG)
gdal_translate -co PROFILE=BASELINE -co COMPRESS=LZW  s2_${AOI}_geotiff_8bit.tif s2_${AOI}_8bit.tif

rm -f ${BAND_R}_geotiff_8bit.tif ${BAND_G}_geotiff_8bit.tif ${BAND_B}_geotiff_8bit.tif s2_${AOI}_geotiff_8bit.tif
Improve this answer
1
  • Thank you very much. – F.H. Jan 27 at 11:24
0

Regarding S2 data as being fully 16-bit will result in "overcompression", as not all the bits within the image are actually used. What you should do is to identify the minimum and maximum values within the raster, and then use those to linearly rescale to a range between 0 and 255.

Basically, that becomes:

rescaled_value=round((pixel_value − min_observed) / (max_observed − min_observed)*255)

And then just save that as 8-bit data.

Improve this answer
7
  • Thanks for answering, but this stretching is not what I'm looking for, as it would make it impossible to compare two images of different epochs, which is necessary to my work. – F.H. Jan 25 at 12:18
  • I need to maintain the relative positions of the DN and be able to see if one scene/band is more energetic than another and stretching the values from 0:255 would hide these crucial information – F.H. Jan 25 at 12:20
  • @F.H. - Then just set min_observed to 1 and max_observed to any sufficiently high value. 65535 is used as the official max value, reserved for saturated pixels. You'll end up with your values highly compressed, once you go to 8 bit, but the relative positions will be maintained. – Mikkel Lydholm Rasmussen Jan 25 at 13:38
  • Thank you, I'll try this approach. I am aware that my data will be heavily compressed, but unfortunately this cannot be avoided for the use of GRASS histogram matching tool, and this is the only tool I have at the moment, till I'll find an alternative that'll allow me to exploit all the information that Sentinel2 imagery contain. – F.H. Jan 25 at 14:41
  • 1
    @F.H. - in most cases, 10000 should be the standard maximum value, but SEN2COR doesn't accurately model many things, including specular reflections, so you may see values significantly above 10k. You could easily cap the values at 10000 and then use that as the maximum value. – Mikkel Lydholm Rasmussen Jan 26 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.