Apply esa_quality_mask

hls_libs/esa_quality_mask/apply_esa_quality_mask.sh

@@ -0,0 +1,47 @@
+# Note:  insert the "esa_quality_mask" code right after "addFmaskSDS" (before "s2trim")
+
+########## What comes before:
+# addFmaskSDS $out_noMASK $fmask20mbin $safexml $granxml $ACCODE $outfile
+
+########## esa_quality_mask
+# Use ESA's image quality mask to filter data.
+# The raster mask is available only for processing baseline 4.0 and above (roughly
+# images produced after January 2022). For earlier processing baseline, nothing to do.
+# 
+# Test if the quality mask exists.
+# Default is that it exists.
+ESA_Quality_Mask=true
+for  band in 01 02 03 04 05 06 07 08 8A 09 10 11 12	
+do 
+	# Convert mask from JP2 to ENVI format for each band. Only record the ENVI B01
+	# filename for later use by the C code, which applies the masks for all bands, 
+	# by inferring mask filename of other bands from B01 mask filename.
+	#
+	# The variable grandir refers to a SAFE directory name.
+	maskjp2=$(ls $grandir/GRANULE/*/QI_DATA/MSK_QUALIT_B${band}.jp2 2>/dev/null)
+	if [ $? -ne 0 ]
+	then
+		# Processing baseline older than 4.0; no quality mask in JP2. 
+		ESA_Quality_Mask=false
+		break
+	fi
+	jp2base=$(basename $maskjp2)
+	maskenvi=$tmpdir/$jp2base.envi
+
+	gdal_translate -of ENVI $maskjp2 $maskenvi
+
+	# Record the filename of the converted B01 quality mask; the filenames of any other band
+	# is inferred from this. 
+	case $band in
+	  01) maskb01=$maskenvi;;
+	esac
+done
+
+# If quality mask exists, apply it. 
+if $ESA_Quality_Mask
+then
+	# Apply mask for all bands; only B01 filename is given as argument.
+	# outfile is the output from the above "addFmaskSDS" and is updated in place.
+	$CODE_BASE_DIR/S2/toTOA/ESA_Quality_Mask/esa_quality_mask $outfile $maskb01
+fi
+

hls_libs/esa_quality_mask/esa_quality_mask.c

@@ -0,0 +1,171 @@
+/******************************************************************************** 
+ * Apply the ESA image quality mask that is available in raster format (JP2) 
+ * (starting with processing baseline 4.0, roughly from Jan 25, 2022 onward) , 
+ * which is much easier to apply than the earlier GML format.  The mask 
+ * indicates whether there is a loss of packets for a pixel in a  spectral 
+ * band, or there is cross-talk, or saturation.
+ *
+ * The details of the mask: 
+ *  https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/data-quality-reports
+ *
+ *  layer   1: lost ancillary packets, 
+ *          2: degraded ancillary packets, 
+ *          3: lost MSI packets, 
+ *          4: degraded MSI packets, 
+ *          5: defective pixels, 
+ *          6: no data, 
+ *          7: partially corrected cross talk, 
+ *          8: saturated pixels. 
+ *  "For each spectral band, the mask is defined at the same spatial resolution as that of the spectral band. 
+ *  In a first step, layer #7 and 8 (partially corrected and saturated) is not activated" 
+ *  (layer 7 is set, contradicting ESA document!)
+ *
+ *  Decided to mask 3 (lost MSI packets) and 4 (degraded MSI packets) only. So far, 5 and 7 are present in 
+ *  B10-B11, because interpolation is used to fill 5 and the cross talk is partially corrected.
+ *  Saturated pixels are still useful for the magnitude of the reflectance; keep it.
+ *
+ * A very late addition, but important.
+ * Sept 20, 2022
+ * Jan 13, 2023
+ * Mar 20, 2023
+ ********************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "hls_projection.h"
+#include "s2r.h"
+#include "util.h"
+
+int main(int argc, char *argv[])
+{
+    /* Commandline parameters */
+    char fname_s2rin[LINELEN];  /* Open hdf for update */
+    char fname_qmB01[LINELEN];   /* B01 quality mask in ENVI format */
+
+    s2r_t s2rin;            
+    int irow, icol, nrow, ncol, k;
+    int npix;
+    int NLAYER = 8; /* Num. of quality layers for each pixel */
+    int psi;        /* pixel size index: 0 for 10m, 1 for 20m, 2 for 60m */
+    int ib, il;     /* band and layer */
+
+    char creationtime[100];
+    char qmname[LINELEN];   /* mask filename for any band */
+    char *base, *cpos;
+    FILE *qmptr;
+    uint8 *mask;
+    char junkchar;      /* A junk character for testing file size */
+    char message[LINELEN];
+
+    char maskSet = 0;       /* Whether the quality mask is set for a band */
+    int8 layer[NLAYER];         /* Which layer of quality mask is set */
+    int ret;
+
+    if (argc != 3) {
+        fprintf(stderr, "%s s2rin  B01_qual_mask.bin \n", argv[0]);
+        exit(1);
+    }
+
+    strcpy(fname_s2rin,  argv[1]);
+    strcpy(fname_qmB01,  argv[2]);
+
+    /* Open the input for update*/
+    strcpy(s2rin.fname, fname_s2rin);
+    ret = open_s2r(&s2rin, DFACC_WRITE);
+    if (ret != 0) {
+        Error("Error in open_s2r");
+        exit(1);
+    }
+
+
+    /* Construct the ESA mask filename from B01 filename */
+    strcpy(qmname, fname_qmB01); 
+    base = strrchr(qmname, '/');
+    if (base == NULL)
+        base = qmname;
+    else
+        base++;
+    /* base is the basename */
+    cpos = strstr(base, "B01");
+    if (cpos == NULL) {
+        fprintf(stderr, "Filename does not contain characters B01: %s\n", fname_qmB01);
+        exit(1);
+    }
+
+    for (ib = 0; ib < S2NBAND; ib++) {
+        strncpy(cpos, S2_SDS_NAME[ib], 3);
+
+        psi = get_pixsz_index(ib);
+        nrow = s2rin.nrow[psi];
+        ncol = s2rin.ncol[psi];
+        npix = nrow*ncol;
+
+        if ((qmptr = fopen(qmname, "r")) == NULL) {
+            sprintf(message, "Can't open for read: %s\n", qmname);
+            Error(message);
+            exit(1);
+        }
+        /* 8 layers for each pixel but some layers are not used; each layer takes 
+         * 1 byte in the ENVI format; originally JP2 must be using a parsimonious 
+         * representation like a bit field, and a bit is turned into a byte during 
+         * conversion by GDAL.
+         */ 
+        if ((mask = (uint8*)calloc(npix*NLAYER, sizeof(uint8))) == NULL) {
+            Error("Can't allocate memory");
+            exit(1);
+        }
+        if (fread(mask, 1, npix*NLAYER, qmptr) != npix*NLAYER ||
+            fread(&junkchar, 1, 1, qmptr) == 1) {
+
+            /* File size should be exact; the extra read attempt is useful guard against
+             * error  */
+            sprintf(message, "File size is wrong: %s\n", qmname);
+            Error(message);
+            exit(1);
+        } 
+
+        maskSet = 0;
+        for (il = 0; il < NLAYER; il++) 
+            layer[il] = -1;
+        for (irow = 0; irow < nrow; irow++) {
+            for (icol = 0; icol < ncol; icol++) {
+                k = irow * ncol + icol;
+                for (il = 0; il < NLAYER; il++) {
+                    if (mask[il*npix+k] == 1) {
+                        if (il == 2 || il == 3)     /* lost or degraded MSI packets; most severe problem */
+                            s2rin.ref[ib][k] = HLS_REFL_FILLVAL;
+
+                        maskSet = 1;
+                        layer[il] = il;     /* Assuming only one layer is set for a given band */
+
+                    }
+
+                }
+            }
+        }
+
+        if (maskSet) { 
+            fprintf(stderr, "maskname = %s\n", qmname); 
+            for (il = 0; il < NLAYER; il++) {
+                if (layer[il] != -1)
+                    fprintf(stderr, "Band %s quality mask is set for layer %d\n", S2_SDS_NAME[ib], layer[il]+1);
+            }
+            fprintf(stderr, "\n");
+        }
+
+        free(mask);
+    }
+
+    /* Processing time */
+    getcurrenttime(creationtime);
+    SDsetattr(s2rin.sd_id, HLSTIME, DFNT_CHAR8, strlen(creationtime), (VOIDP)creationtime);
+
+    if (close_s2r(&s2rin) != 0) {
+        Error("Error in closing output");
+        exit(1);
+    }
+
+    return 0;
+}

hls_libs/esa_quality_mask/makefile

@@ -0,0 +1,49 @@
+# Note: before run make, do in the shell:
+# source $HOME/code.toa/share/makeenv/makeenv.sh local
+
+TGT = esa_quality_mask
+OBJ = 	esa_quality_mask.o \
+	hls_projection.o \
+	s2r.o \
+	util.o \
+	hdfutility.o 
+
+$(TGT): $(OBJ)
+	$(CC) $(CFLAGS) -o $(TGT) $(OBJ) -L$(GCTPLIB) -L$(HDFLIB) -L$(ZLIB) -L$(SZLIB) -L$(JPGLIB) -L$(PROJLIB)  $(GCTPLINK) $(HDFLINK) 
+
+esa_quality_mask.o: esa_quality_mask.c 
+	$(CC) $(CFLAGS) -c esa_quality_mask.c -I$(HDFINC) -I$(GCTPINC)
+
+hls_projection.o: hls_projection.c
+	$(CC) $(CFLAGS) -c hls_projection.c -I$(HDFINC) -I$(GCTPINC) 
+
+s2r.o: s2r.c
+	$(CC) $(CFLAGS) -c s2r.c -I$(HDFINC) -I$(GCTPINC)
+
+util.o: util.c
+	$(CC) $(CFLAGS) -c  util.c -I$(HDFINC)
+
+hdfutility.o: hdfutility.c
+	$(CC) $(CFLAGS) -c  hdfutility.c -I$(HDFINC)
+
+link:
+	ln -sf ../../common/s2r.h
+	ln -sf ../../common/s2r.c
+	ln -sf ../../common/s2def.h
+	ln -sf ../../common/s2at30m.h		# included by hls_hdfeos.h 	
+	ln -sf ../../../L8/common/lsat.h
+	ln -sf ../../../share/common/cubic_conv.h 	# included by a2at30m.h
+	ln -sf ../../../share/common/util.h
+	ln -sf ../../../share/common/util.c
+	ln -sf ../../../share/common/hdfutility.h
+	ln -sf ../../../share/common/hdfutility.c
+	ln -sf ../../../share/common/hls_projection.h
+	ln -sf ../../../share/common/hls_projection.c
+	ln -sf ../../../share/common/hls_hdfeos.h
+	ln -sf ../../../share/common/hls_hdfeos.c
+	ln -sf ../../../share/common/fillval.h
+	ln -sf ../../../share/common/hls_commondef.h
+
+clean:
+	rm -f *.o
+