New code for S2 twin granule consolidation

Dockerfile

@@ -111,9 +111,9 @@ RUN pip3 install --upgrade awscli
 RUN pip3 install click==7.1.2
 RUN pip3 install rio-cogeo==1.1.10 --no-binary rasterio --user
 
-RUN pip3 install git+https://github.com/NASA-IMPACT/hls-thumbnails@v1.1
+RUN pip3 install git+https://github.com/NASA-IMPACT/hls-thumbnails@v1.3
 
-RUN pip3 install git+https://github.com/NASA-IMPACT/hls-metadata@v2.2
+RUN pip3 install git+https://github.com/NASA-IMPACT/hls-metadata@v2.4
 
 RUN pip3 install git+https://github.com/NASA-IMPACT/[email protected]
 

hls_libs/consolidate/consolidate.c

@@ -20,6 +20,11 @@
  *
  * Apr 1, 2020: Note that the mean sun/view angles are from twin A only; later in 
  *   derive_s2nbar these angles will be recomputed and so set properly. 
+ *
+ * Mar 18, 2024:
+ *	 The silly max NDVI rule in selecting between two twin granules is eliminated
+ *	 because the leading/trailing edges of the twin granules have been trimmed where
+ *	 some pixels are possibly corrupted. With trimming, data from either granule is good. 
  ********************************************************************************/
 
 #include "s2r.h"
@@ -31,9 +36,6 @@ int copy_metadata_AB(s2r_t *s2rA, s2r_t *s2rB, s2r_t *s2rO);
 /* Keep the pixel whose 10m row/col are given */
 int copypix(s2r_t *from, int irow60m, int icol60m, s2r_t *to);
 
-/* Compute NDVI for the 60m pixel location */
-double ndvi(s2r_t *s2, int irow60m, int icol60m);
-
 int main(int argc, char * argv[])
 {
 	/* Commandline parameters */
@@ -99,12 +101,8 @@ int main(int argc, char * argv[])
 			if (s2rA.ref[ubidx][k60m] != HLS_REFL_FILLVAL ) 
 				copypix(&s2rA, irow60m, icol60m, &s2rO);
 
-			if (s2rB.ref[ubidx][k60m] != HLS_REFL_FILLVAL) {
-				if (s2rO.ref[ubidx][k60m] == HLS_REFL_FILLVAL) 
-					copypix(&s2rB, irow60m, icol60m, &s2rO);
-				else if (ndvi(&s2rB, irow60m, icol60m) > ndvi(&s2rO, irow60m,icol60m)) 
+			if (s2rB.ref[ubidx][k60m] != HLS_REFL_FILLVAL && s2rO.ref[ubidx][k60m] == HLS_REFL_FILLVAL) 
 					copypix(&s2rB, irow60m, icol60m, &s2rO);
-			}
 		}
 	}
 
@@ -263,34 +261,3 @@ int copypix(s2r_t *from, int irow60m, int icol60m, s2r_t *to)
 	return 0;
 }
 
-/* Compute NDVI for the 60m pixel location */
-double ndvi(s2r_t *s2, int irow60m, int icol60m)
-{
-	int irow, icol;
-	int ncol, k;
-	int rowbeg, rowend, colbeg, colend;
-	double red, nir;
-	int bs = 6; 	/* 6 10m pixels nesting within a 60m pixel in one dimension */
-
-	ncol = s2->ncol[0]; 	/* 10m pixels */
-
-	rowbeg = irow60m * bs;
-	rowend = rowbeg + bs - 1;
-	colbeg = icol60m * bs;
-	colend = colbeg + bs - 1;
-
-	red = nir = 0;
-	for (irow = rowbeg; irow <= rowend; irow++) {
-		for (icol = colbeg; icol <= colend; icol++) {
-			k = irow * ncol + icol;
-			if (s2->ref[3][k] != HLS_REFL_FILLVAL && s2->ref[7][k] != HLS_REFL_FILLVAL) {
-      				red += s2->ref[3][k];
-				nir += s2->ref[7][k];
-			}
-		}
-	}
-	if (red != 0 && nir != 0)
-		return (nir-red)/(nir+red);
-	else
-		return HLS_REFL_FILLVAL;
-}

hls_libs/trim/s2trim.c

@@ -1,9 +1,20 @@
 /******************************************************************************** 
- * Trim the image edge so that in the output a location either has data in all 
- * spectral bands or has no data in any spectral band. 
+ * 1) The images of different spectral bands do not start from and end at the 
+ * same ground location. Trim the image edge so that in the output a location 
+ * either has data in all spectral bands or has no data in any spectral band. 
+ * The trimming has effect in the interior of the images too. If atmospheric
+ * correction makes the data invalid at a pixel in one spectral band, the pixel
+ * will have no data in all spectral bands after trimming.
+ * 
+ * 2) The leading or trailing edge of a twin granule can have corrupted data that
+ * ESA does not flag. So trim the edge in all spectral bands by a length of n
+ * 60-m pixels. N = 2 for now.
+ * The side effect is that the leading/trailing edge of the data take will be 
+ * trimmed too, but the data loss there has no significance.
+ *
+ * First implemented on Apr 8, 2020, and added trailing/leading edge trimming 
+ * on Mar 4, 2024
  *
- * A very late addition.
- * Apr 8, 2020
  ********************************************************************************/
 
 #include <stdio.h>
@@ -27,6 +38,12 @@ int main(int argc, char *argv[])
 	char missing;
 	int bs;
 
+	/*** A few variables used for leading/trailing edge trimming */
+	char isFill;
+	/* Number of pixels to be trimmed in the 60-m bands at the leading/trailing edge*/
+	int N60m = 2;	
+	int startRow, endRow; /* Defining the rows to be trimmed */
+
 	char creationtime[100];
 	int ret;
 
@@ -45,12 +62,81 @@ int main(int argc, char *argv[])
 		exit(1);
 	}
 
-	/* Use the 60m aerosol band to guide the trimming. For a 60m by 60m area if
-	 * there is no measurement in any spectral band, the entire 60m by 60m 
-	 * area will filled with nodata.
+	/* Use the 60m aerosol band to guide the trimming. If there is nodata value in any 
+	 * spectral band for the 60m x 60m pixel, the entire 60m x 60m area in all spectral 
+	 * bands area will be filled with nodata.
 	 */
 	nrow60m = s2rin.nrow[2];
 	ncol60m = s2rin.ncol[2];
+
+	/* Leading/trailing edge trimming in the coastal/aerosol band by N pixels, and this
+	 * will have an effect on the subsquent trimming in all bands. Added March 2024.
+	 */
+	for (icol60m = 0; icol60m < ncol60m; icol60m++) {
+
+		/*** First check if the top of each column has fill value that indicates the edge */
+		isFill = 0;	/* Fillvalue or not */
+		startRow = -1;	/* Looking from the top down. The row that starts to have non-fill 
+						 * value and the previous row has fill value */
+		for (irow60m = 0; irow60m < nrow60m; irow60m++) {  
+			k60m = irow60m * ncol60m + icol60m;
+			if (s2rin.ref[0][k60m] == HLS_REFL_FILLVAL) {
+				isFill = 1;
+				continue;
+			}
+			else {
+				if (isFill == 1)	/* The row above is fill; edge detected */
+					startRow = irow60m;
+
+				break;	/* Break whenever nonfill is found */
+			}
+		}
+
+		/* Set pixels in this column from startRow to endRow to fill value */
+		if (startRow != -1) {
+			endRow = startRow + N60m-1;
+			if (endRow >= nrow60m)
+				endRow = nrow60m-1;	/* the last row of the image in the 60m band*/
+
+			for (irow60m = startRow; irow60m <= endRow; irow60m++) {  /* including endRow */
+				k60m = irow60m * ncol60m + icol60m;
+				s2rin.ref[0][k60m] = HLS_REFL_FILLVAL;
+			}
+		}
+
+		/*** Then check if the bottom of each column has fill value */
+		isFill = 0;	/* Fillvalue or not */
+		endRow = -1;	/* Looking from the bottom up. The row that starts to have non-fill 
+						 * value and the previous row is fill value */
+		for (irow60m = nrow60m-1; irow60m >= 0; irow60m--) {  
+			k60m = irow60m * ncol60m + icol60m;
+			if (s2rin.ref[0][k60m] == HLS_REFL_FILLVAL) {
+				isFill = 1;
+				continue;
+			}
+			else {
+				if (isFill == 1)	/* The row below is fill; edge detected */
+					endRow = irow60m;
+
+				break;	/* Break whenever nonfill is found */
+			}
+		}
+
+		/* Set pixels in this column from startRow to endRow to fill value */
+		if (endRow != -1) {
+			startRow = endRow - (N60m-1);
+			if (startRow <= 0)
+				startRow = 0;
+
+			for (irow60m = startRow; irow60m <= endRow; irow60m++) {  /* including endRow */
+				k60m = irow60m * ncol60m + icol60m;
+				s2rin.ref[0][k60m] = HLS_REFL_FILLVAL;
+			}
+		}
+	}	/* End of leading/trailing edge trimming */
+
+
+
 	for (irow60m = 0; irow60m < nrow60m; irow60m++) {
 		for (icol60m = 0; icol60m < ncol60m; icol60m++) {
 			k60m = irow60m * ncol60m + icol60m;