From f0755f02fe8be9bc803d2bf52beb46dd2a5576da Mon Sep 17 00:00:00 2001 From: junchangju Date: Wed, 20 Mar 2024 12:05:45 -0400 Subject: [PATCH 1/3] New code for S2 twin granule consolidation --- hls_libs/consolidate/consolidate.c | 45 ++----------- hls_libs/trim/s2trim.c | 100 +++++++++++++++++++++++++++-- 2 files changed, 99 insertions(+), 46 deletions(-) diff --git a/hls_libs/consolidate/consolidate.c b/hls_libs/consolidate/consolidate.c index 8a24309..a80e744 100755 --- a/hls_libs/consolidate/consolidate.c +++ b/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; -} diff --git a/hls_libs/trim/s2trim.c b/hls_libs/trim/s2trim.c index 3dfa365..cc5abc4 100755 --- a/hls_libs/trim/s2trim.c +++ b/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 @@ -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; From 6d938fb9bcfddbe2cf264af8ee3093f5e4b997f0 Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Thu, 2 May 2024 13:44:51 -0600 Subject: [PATCH 2/3] Update hls-thumbnails version to avoid broken build. --- Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index 78e5e0a..3799989 100644 --- a/Dockerfile +++ b/Dockerfile @@ -111,7 +111,7 @@ 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 From 2590ddb6061583653ba2c39074b9bf954a3452dc Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Fri, 3 May 2024 20:01:48 -0600 Subject: [PATCH 3/3] Update hls-metadata version to avoid broken build. --- Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index 3799989..c8593c5 100644 --- a/Dockerfile +++ b/Dockerfile @@ -113,7 +113,7 @@ RUN pip3 install rio-cogeo==1.1.10 --no-binary rasterio --user 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/hls-manifest@v2.0