update after paper 2 is accepted

fsn1995 · Dec 30, 2023 · 0e6ed4a · 0e6ed4a
1 parent 2e9a26e
commit 0e6ed4a
Show file tree

Hide file tree

Showing 16 changed files with 473 additions and 50 deletions.
diff --git a/QGIS/albedoQGIS.py b/QGIS/albedoQGIS.py
@@ -247,7 +247,7 @@ def prepS2(img):
             .filter(colFilter) \
             .map(prepEtm) \
             .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])
-s2Col = ee.ImageCollection('COPERNICUS/S2_SR') \
+s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \
             .filter(s2colFilter) \
             .map(prepS2) \
             .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])

diff --git a/QGIS/arcticdemQGIS.py b/QGIS/arcticdemQGIS.py
@@ -20,4 +20,27 @@
 
 #%% GLIMS
 glims = ee.FeatureCollection('GLIMS/current')
-Map.addLayer(arcticDEM.select('elevation').clip(glims), visPara, 'arctic dem glims')
+Map.addLayer(arcticDEM.select('elevation').clip(glims), visPara, 'arctic dem glims')
+
+
+#%% make colorbar
+# ref https://matplotlib.org/stable/tutorials/colors/colorbar_only.html
+# import matplotlib.pyplot as plt
+# import matplotlib as mpl
+# from matplotlib.colors import ListedColormap#, LinearSegmentedColormap
+
+
+# fig, ax = plt.subplots(figsize=(6, 1))
+# fig.subplots_adjust(bottom=0.5)
+
+# blue_fluorite = ['#291b32', '#2a1b34', '#2b1b34', '#2d1c36', '#2f1c38', '#301c39', '#301d3a', '#321d3b', '#331d3d', '#351d3f', '#351e40', '#371e41', '#381e43', '#3a1e45', '#3b1f45', '#3c1f46', '#3e1f48', '#3f1f4a', '#401f4c', '#42204d', '#43204e', '#44204f', '#462051', '#472052', '#482054', '#4a2056', '#4a2157', '#4c2158', '#4e215a', '#4f215b', '#50215d', '#52215e', '#532160', '#552162', '#552263', '#562264', '#582265', '#592267', '#5b2268', '#5c226b', '#5e226c', '#5f226e', '#60226f', '#622271', '#632272', '#642274', '#662276', '#672277', '#692278', '#6a227a', '#6c227b', '#6e227d', '#6e237e', '#6f247f', '#702480', '#712581', '#722681', '#732683', '#742783', '#752884', '#762985', '#772987', '#792a87', '#792b88', '#7a2c89', '#7b2c8a', '#7c2d8a', '#7d2d8c', '#7e2e8d', '#7f2f8d', '#80308e', '#813190', '#823191', '#833292', '#843292', '#863393', '#863494', '#873595', '#893596', '#8a3697', '#8b3798', '#8b3899', '#8c389a', '#8e399b', '#8e3a9c', '#8f3b9c', '#8f3d9d', '#8f3e9e', '#903f9e', '#90419e', '#90439f', '#9044a0', '#9046a0', '#9047a1', '#9049a1', '#914aa2', '#914ca2', '#914ca3', '#914ea3', '#9150a4', '#9151a5', '#9153a5', '#9154a6', '#9156a6', '#9157a7', '#9258a7', '#9259a8', '#925aa8', '#925ba9', '#925da9', '#925faa', '#9260ab', '#9260ab', '#9263ac', '#9264ac', '#9265ad', '#9266ae', '#9268ae', '#9269ae', '#926aaf', '#926bb0', '#926cb0', '#926eb1', '#926fb1', '#9270b2', '#9271b2', '#9273b3', '#9274b3', '#9275b4', '#9277b5', '#9277b5', '#9278b6', '#927ab6', '#927bb7', '#927cb7', '#927eb8', '#927fb8', '#9280b9', '#9281ba', '#9282ba', '#9284bb', '#9285bb', '#9285bc', '#9187bc', '#9188bd', '#918abd', '#918bbe', '#918cbf', '#918dbf', '#918ec0', '#918fc0', '#9191c1', '#9092c2', '#9094c2', '#9094c2', '#9095c3', '#9096c3', '#8f99c4', '#8f9ac5', '#8f9ac5', '#8f9bc6', '#8f9cc6', '#8f9dc7', '#8e9fc8', '#8ea0c8', '#8ea2c9', '#8ea3c9', '#8da5ca', '#8da5ca', '#8da6cb', '#8da7cb', '#8ca9cc', '#8caacc', '#8caccd', '#8bacce', '#8badce', '#8baecf', '#8ab0d0', '#8ab2d0', '#8ab2d1', '#8ab4d1', '#89b4d1', '#89b5d2', '#89b7d2', '#88b8d3', '#88bad4', '#87bad4', '#87bbd5', '#86bdd6', '#86bed6', '#86c0d7', '#85c0d7', '#85c1d8', '#84c3d8', '#84c4d9', '#83c5d9', '#83c6da', '#82c8da', '#82c8db', '#81cadc', '#81cbdc', '#80ccdd', '#81cddd', '#84cfdd', '#85cfdd', '#87d0dd', '#8ad0de', '#8dd1de', '#8fd2de', '#90d2de', '#92d4de', '#95d5de', '#97d5de', '#98d6de', '#9bd7de', '#9dd7df', '#a0d8df', '#a1d9df', '#a2dadf', '#a5dadf', '#a7dbdf', '#aadcdf', '#abdddf', '#acdde0', '#afdfe0', '#b1dfe0', '#b3e0e0', '#b4e1e0', '#b7e2e0', '#bae2e1', '#bae3e1', '#bee3e2', '#c0e4e3', '#c1e5e3', '#c4e6e3', '#c6e6e4', '#c8e7e4', '#cbe7e5', '#cde8e5', '#cee9e6', '#d2e9e7', '#d3eae7', '#d5eae7', '#d8ebe8', '#d9ece8', '#dcece9', '#deedea', '#dfeeea', '#e2eeea', '#e5efeb', '#e6f0eb', '#e9f0ec', '#ebf1ed', '#ecf2ed', '#eff3ee', '#f1f3ee']
+# cmp = ListedColormap(blue_fluorite)
+# cmap = mpl.cm.cool
+# norm = mpl.colors.Normalize(vmin=0, vmax=1)
+
+# cb1 = mpl.colorbar.ColorbarBase(ax, cmap=cmp,
+#                                 norm=norm,
+#                                 orientation='horizontal')
+# cb1.set_label('albedo')
+# # fig.show()
+# fig.savefig("colorbar.svg", dpi=300, transparent=True)
diff --git a/README.md b/README.md
@@ -5,8 +5,8 @@
 
 This is a github repo for remote sensing of ice/snow albedo using Google Earth Engine.
 It is part of two studies:
-1. <div class="csl-entry">Feng, S., Cook, J. M., Anesio, A. M., Benning, L. G., &#38; Tranter, M. (2023). Long time series (1984–2020) of albedo variations on the Greenland ice sheet from harmonized Landsat and Sentinel 2 imagery. <i>Journal of Glaciology</i>, 1–16. https://doi.org/10.1017/jog.2023.11</div>
-2. Remote sensing of ice albedo using harmonized Landsat and Sentinel 2 datasets: Validation. The manuscript is submitted for peer review. 
+1. Feng, S., Cook, J. M., Anesio, A. M., Benning, L. G., & Tranter, M. (2023). Long time series (1984–2020) of albedo variations on the Greenland ice sheet from harmonized Landsat and Sentinel 2 imagery. Journal of Glaciology, 69(277), 1225–1240. [https://doi.org/10.1017/jog.2023.11](https://doi.org/10.1017/jog.2023.11).
+2. Feng, S., Cook, J. M., Onuma, Y., Naegeli, K., Tan, W., Anesio, A. M., Benning, L. G., & Tranter, M. (2023). Remote sensing of ice albedo using harmonized Landsat and Sentinel 2 datasets: validation. International Journal of Remote Sensing, 00(00), 1–29. [https://doi.org/10.1080/01431161.2023.2291000](https://doi.org/10.1080/01431161.2023.2291000).
 
 Paper 1 focuses on the development of data harmonization and narrow to broad band conversion algorithm optimized for the ablation zone in the Greenland Ice Sheet. The harmonized satellite albedo (HSA) is further validated globally using AWS albedo measurements in paper 2. It also provides a complete workflow of validating remote sensing derived product with groundtruth measurements. 
 
@@ -23,6 +23,12 @@ This repository relates to the harmonization of Landsat and Sentinel-2 data and
 Here is a simple tutorial that demonstrates the importance of harmonizing Landsat 4-7, 9 and Sentinel 2 to Landsat 8 time series of datasets.
 It will display the charts of the harmonized satellite albedo (All Observations) and original albedo (All Observations Original).
 The linear trendline will be plotted on a separate chart. 
+It can be used to determine whether data harmonziation is necessary or not, and to extract the time series of albedo at a specific location.
+
+-[script/hsaImgFinder.js](script\hsaImgFinder.js)
+
+This is a simple tutorial that demonstrates how to find the harmonized satellite albedo image for a specific date and location.
+It will display the rgb and albedo image. It supports exporting the image to Google Drive as well.
 
 <img src="media\harmonization.png" alt="harmonization" height=800/>
 
@@ -183,7 +189,8 @@ This project is organized into three main directories. `geemap` contains the jav
 |   |   darkzone.ipynb
 |   |   darkzonePoint.ipynb
 |   |   geopyfsn.py
-|   |   harmonizationTutorial.js        
+|   |   harmonizationTutorial.js       
+|   |   hsaImgFinder.js 
 |   |   L7ToL8.js
 |   |   L9ToL8.js
 |   |   Promice vs satellite modis.ipynb

diff --git a/geeapp/albedoInspectorAPP.js b/geeapp/albedoInspectorAPP.js
@@ -63,8 +63,8 @@ function renameEtm(img) {
 // Function to get and rename bands of interest from Sentinel 2.
 function renameS2(img) {
   return img.select(
-    ['B2',   'B3',    'B4',  'B8',  'QA60', 'SCL'],
-    ['Blue', 'Green', 'Red', 'NIR', 'QA60', 'SCL']
+    ['B2',   'B3',    'B4',  'B8',  'QA60', 'SCL', QA_BAND],
+    ['Blue', 'Green', 'Red', 'NIR', 'QA60', 'SCL', QA_BAND]
     //['B2',     'B3',      'B4',    'B8',    'B11',     'B12',     'QA60', 'SCL'],
     //['BlueS2', 'GreenS2', 'RedS2', 'NIRS2', 'SWIR1S2', 'SWIR2S2', 'QA60', 'SCL']
   );
@@ -163,21 +163,44 @@ function maskL457sr(image) {
  * Function to mask clouds using the Sentinel-2 QA band
  * @param {ee.Image} image Sentinel-2 image
  * @return {ee.Image} cloud masked Sentinel-2 image
+ * archived after updating to Cloud Score+
  */
-function maskS2sr(image) {
-  var qa = image.select('QA60');
+// function maskS2sr(image) {
+//   var qa = image.select('QA60');
+
+//   // Bits 10 and 11 are clouds and cirrus, respectively.
+//   var cloudBitMask = 1 << 10;
+//   var cirrusBitMask = 1 << 11;
+//   // 1 is saturated or defective pixel
+//   var not_saturated = image.select('SCL').neq(1);
+//   // Both flags should be set to zero, indicating clear conditions.
+//   var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
+//       .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
+
+//   // return image.updateMask(mask).updateMask(not_saturated);
+//   return image.updateMask(mask).updateMask(not_saturated).divide(10000);
+// }
+
+/**
+ * Function to mask clouds using the Cloud Score+
+ */
+// Cloud Score+ image collection. Note Cloud Score+ is produced from Sentinel-2
+// Level 1C data and can be applied to either L1C or L2A collections.
+var csPlus = ee.ImageCollection('GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED');
 
-  // Bits 10 and 11 are clouds and cirrus, respectively.
-  var cloudBitMask = 1 << 10;
-  var cirrusBitMask = 1 << 11;
+// Use 'cs' or 'cs_cdf', depending on your use case; see docs for guidance.
+var QA_BAND = 'cs'; // I find'cs' is better than 'cs_cdf' because it is more robust but may mask out more clear pixels.
+
+// The threshold for masking; values between 0.50 and 0.65 generally work well.
+// Higher values will remove thin clouds, haze & cirrus shadows.
+var CLEAR_THRESHOLD = 0.65;
+
+function maskS2sr(image) {
   // 1 is saturated or defective pixel
   var not_saturated = image.select('SCL').neq(1);
-  // Both flags should be set to zero, indicating clear conditions.
-  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
-      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
-
-  // return image.updateMask(mask).updateMask(not_saturated);
-  return image.updateMask(mask).updateMask(not_saturated).divide(10000);
+  return image.updateMask(image.select(QA_BAND).gte(CLEAR_THRESHOLD))
+              .updateMask(not_saturated)
+              .divide(10000);
 }
 
 // narrow to broadband conversion
@@ -342,7 +365,8 @@ var generateChart = function (coords) {
     ee.Filter.bounds(point),
     ee.Filter.date(date_start, date_end),
     // ee.Filter.calendarRange(6, 7, 'month'),
-    ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 50)
+    // ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 50),
+    ee.Filter.lt('MEAN_SOLAR_ZENITH_ANGLE', 76)
   );
 
   var oli2Col = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2') 
@@ -369,6 +393,7 @@ var generateChart = function (coords) {
               .select(['visnirAlbedo']); //# .select(['totalAlbedo']) or  .select(['visnirAlbedo'])
   var s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') 
               .filter(s2colFilter) 
+              .linkCollection(csPlus, [QA_BAND])
               .map(prepS2)
               .select(['visnirAlbedo']); //# .select(['totalAlbedo']) or  .select(['visnirAlbedo'])
 
@@ -592,7 +617,8 @@ var loadComposite = function() {
       ee.Filter.bounds(aoi),
       ee.Filter.date(startDate, endDate),
       // ee.Filter.calendarRange(6, 7, 'month'),
-      ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 50)
+      // ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 50),
+      ee.Filter.lt('MEAN_SOLAR_ZENITH_ANGLE', 76)
     );
 
     var oli2Col = ee.ImageCollection('LANDSAT/LC09/C02/T1_L2') 
@@ -618,6 +644,7 @@ var loadComposite = function() {
                 .map(prepEtm);
                 // .select(['visnirAlbedo']); //# .select(['totalAlbedo']) or  .select(['visnirAlbedo'])
     var s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') 
+                .linkCollection(csPlus, [QA_BAND])
                 .filter(s2colFilter) 
                 .map(prepS2);
                 // .select(['visnirAlbedo']); //# .select(['totalAlbedo']) or  .select(['visnirAlbedo'])
@@ -683,9 +710,10 @@ var logoPanel = ui.Panel(thumb, 'flow', {width: '120px'});
 inspectorPanel.widgets().set(6, logoPanel);
 
 var logoIntro = ui.Panel([
-  ui.Label("The Deep Purple project receives funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 856416. This study is currently under review."),
+  ui.Label("The Deep Purple project receives funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 856416."),
   ui.Label("https://www.deeppurple-ercsyg.eu/home", {}, "https://www.deeppurple-ercsyg.eu/home"),
   ui.Label("https://github.com/fsn1995/Remote-Sensing-of-Albedo", {}, "https://github.com/fsn1995/Remote-Sensing-of-Albedo"),
-  ui.Label("Feng, S., Cook, J. M., Anesio, A. M., Benning, L. G. and Tranter, M. (2023) “Long time series (1984–2020) of albedo variations on the Greenland ice sheet from harmonized Landsat and Sentinel 2 imagery,” Journal of Glaciology. Cambridge University Press, pp. 1–16. doi: 10.1017/jog.2023.11.", {}, "https://doi.org/10.1017/jog.2023.11")
+  ui.Label("Feng, S., Cook, J. M., Anesio, A. M., Benning, L. G. and Tranter, M. (2023) “Long time series (1984–2020) of albedo variations on the Greenland ice sheet from harmonized Landsat and Sentinel 2 imagery,” Journal of Glaciology. Cambridge University Press, 69(277), pp. 1225–1240. doi: 10.1017/jog.2023.11.", {}, "https://doi.org/10.1017/jog.2023.11"),
+  ui.Label("Feng, S., Cook, J. M., Onuma, Y., Naegeli, K., Tan, W., Anesio, A. M., Benning, L. G., & Tranter, M. (2023). Remote sensing of ice albedo using harmonized Landsat and Sentinel 2 datasets: validation. International Journal of Remote Sensing, 00(00), 1–29. https://doi.org/10.1080/01431161.2023.2291000", {}, "https://doi.org/10.1080/01431161.2023.2291000")
 ]);
 inspectorPanel.widgets().set(7, logoIntro);
diff --git a/script/L7ToL8.js b/script/L7ToL8.js
@@ -162,7 +162,7 @@ var etmCol = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
 // var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
 //             .filter(colFilter)
 //             .map(prepEtm);
-// var s2Col = ee.ImageCollection('COPERNICUS/S2_SR')
+// var s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
 //             .filter(s2colFilter)
 //             .map(prepS2);
 

diff --git a/script/L9ToL8.js b/script/L9ToL8.js
@@ -136,7 +136,7 @@ var oli2Col = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
 // var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
 //             .filter(colFilter)
 //             .map(prepEtm);
-// var s2Col = ee.ImageCollection('COPERNICUS/S2_SR')
+// var s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
 //             .filter(s2colFilter)
 //             .map(prepS2);
 

diff --git a/script/S2TOL8.js b/script/S2TOL8.js
@@ -184,7 +184,7 @@ var oliCol = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
              .filter(colFilter)
              .map(prepOli);
 
-var s2Col = ee.ImageCollection('COPERNICUS/S2_SR')
+var s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
             .filter(s2colFilter)
             .map(prepS2);
 

diff --git a/script/darkzone.ipynb b/script/darkzone.ipynb
@@ -430,7 +430,7 @@
     "            .filter(colFilter) \\\n",
     "            .map(prepEtm) \\\n",
     "            .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])\n",
-    "s2Col = ee.ImageCollection('COPERNICUS/S2_SR') \\\n",
+    "s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \\\n",
     "            .filter(s2colFilter) \\\n",
     "            .map(prepS2) \\\n",
     "            .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])\n",
@@ -726,7 +726,7 @@
     "            .filter(colFilter) \\\n",
     "            .map(prepEtm) \\\n",
     "            .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])\n",
-    "s2Col = ee.ImageCollection('COPERNICUS/S2_SR') \\\n",
+    "s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \\\n",
     "            .filter(s2colFilter) \\\n",
     "            .map(prepS2) \\\n",
     "            .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])\n",

diff --git a/script/darkzonePoint.ipynb b/script/darkzonePoint.ipynb
@@ -305,7 +305,7 @@
     "            .filter(colFilter) \\\n",
     "            .map(prepEtm) \\\n",
     "            .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])\n",
-    "s2Col = ee.ImageCollection('COPERNICUS/S2_SR') \\\n",
+    "s2Col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \\\n",
     "            .filter(s2colFilter) \\\n",
     "            .map(prepS2) \\\n",
     "            .select(['visnirAlbedo']) # .select(['totalAlbedo']) or  .select(['viznirAlbedo'])\n",