Merge pull request #5 from emit-sds/develop

Merge develop into main for v1.3.0
emit-sds · Nov 7, 2022 · d3d6e3a · d3d6e3a
2 parents 80acdd5 + 2953192
commit d3d6e3a
Show file tree

Hide file tree

Showing 6 changed files with 264 additions and 30 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,13 +4,24 @@ All notable changes to this project will be documented in this file. Dates are d
 
 Generated by [`auto-changelog`](https://github.com/CookPete/auto-changelog).
 
+#### [v1.3.0](https://github.com/emit-sds/emit-sds-l2a/compare/v1.2.0...v1.3.0)
+
+> 3 November 2022
+
+- Analytical atm updates [`#4`](https://github.com/emit-sds/emit-sds-l2a/pull/4)
+- Rerun updates [`#3`](https://github.com/emit-sds/emit-sds-l2a/pull/3)
+- updated covariance for oxygen a help [`272d852`](https://github.com/emit-sds/emit-sds-l2a/commit/272d852d89bac2dab5adc597deebf1d16b595d1f)
+- updated surface priors [`692e676`](https://github.com/emit-sds/emit-sds-l2a/commit/692e676d3c41886904e4f697e614b0efeab23049)
+- add new surface model [`ca436d2`](https://github.com/emit-sds/emit-sds-l2a/commit/ca436d2b1fde3617e9e0d58ccfe57a004a8318a2)
+
 #### [v1.2.0](https://github.com/emit-sds/emit-sds-l2a/compare/v1.1.0...v1.2.0)
 
-> 18 July 2022
+> 22 July 2022
 
+- Merge develop into main for v1.2.0 [`#2`](https://github.com/emit-sds/emit-sds-l2a/pull/2)
 - renamed main directory, and new config [`b236d96`](https://github.com/emit-sds/emit-sds-l2a/commit/b236d9685edf336a983c23ff1762c86b0af50a67)
 - renamed, and removed some cruft [`5e894b5`](https://github.com/emit-sds/emit-sds-l2a/commit/5e894b5d03bc6580eec5c68caf4ea59774f01451)
-- include alternate keystrings for aerosols (facilitates emulator) [`2337c93`](https://github.com/emit-sds/emit-sds-l2a/commit/2337c93c2f4a0f900cdd658950abc68908bfc39f)
+- Update change log [`1b3d9cc`](https://github.com/emit-sds/emit-sds-l2a/commit/1b3d9ccf6824cba48ce84705490188f1e441004d)
 
 #### [v1.1.0](https://github.com/emit-sds/emit-sds-l2a/compare/v1.0.0...v1.1.0)
 

diff --git a/make_emit_masks.py b/make_emit_masks.py
@@ -44,7 +44,7 @@ def haversine_distance(lon1, lat1, lon2, lat2, radius=6335439):
 
 
 @ray.remote
-def build_line_masks(start_line: int, stop_line: int, rdnfile: str, locfile: str, lblfile: str, state: np.array, dt: datetime, h2o_band: np.array, aod_bands: np.array, pixel_size: float, outfile: str, wl: np.array, irr: np.array):
+def build_line_masks(start_line: int, stop_line: int, rdnfile: str, locfile: str, atmfile: str, dt: datetime, h2o_band: np.array, aod_bands: np.array, pixel_size: float, outfile: str, wl: np.array, irr: np.array):
     # determine glint bands having negligible water reflectance
     BLUE = np.logical_and(wl > 440, wl < 460)
     NIR = np.logical_and(wl > 950, wl < 1000)
@@ -61,15 +61,14 @@ def build_line_masks(start_line: int, stop_line: int, rdnfile: str, locfile: str
 
     rdn_ds = envi.open(envi_header(rdnfile)).open_memmap(interleave='bil')
     loc_ds = envi.open(envi_header(locfile)).open_memmap(interleave='bil')
-    lbl_ds = envi.open(envi_header(lblfile)).open_memmap(interleave='bil')
+    atm_ds = envi.open(envi_header(atmfile)).open_memmap(interleave='bil')
 
     return_mask = np.zeros((stop_line - start_line, 8, rdn_ds.shape[2]))
     for line in range(start_line, stop_line):
         print(f'{line} / {stop_line - start_line}')
         loc = loc_ds[line,...].copy().astype(np.float32).T
         rdn = rdn_ds[line,...].copy().astype(np.float32).T
-        lbl = lbl_ds[line,...].copy().astype(np.float32).T
-        x = np.zeros((rdn.shape[0], state.shape[1]))
+        atm = atm_ds[line,...].copy().astype(np.float32).T
 
         elevation_m = loc[:, 2]
         latitude = loc[:, 1]
@@ -79,8 +78,8 @@ def build_line_masks(start_line: int, stop_line: int, rdnfile: str, locfile: str
 
         rho = (((rdn * np.pi) / (irr.T)).T / np.cos(zen)).T
 
-        rho[rho[:, 0] < -9990, :] = -9999.0
-        bad = (latitude < -9990).T
+        rho[rho[:, 0] <= -9990, :] = -9999.0
+        bad = (latitude <= -9990).T
 
         # Cloud threshold from Sandford et al.
         total = np.array(rho[:, b450] > 0.28, dtype=int) + \
@@ -100,23 +99,17 @@ def build_line_masks(start_line: int, stop_line: int, rdnfile: str, locfile: str
         # Threshold spacecraft parts using their lack of an O2 A Band
         mask[3, :] = np.array(rho[:, b762]/rho[:, b780] > 0.8, dtype=int)
 
-        for i, j in enumerate(lbl[:, 0]):
-            if j <= 0:
-                x[i, :] = -9999.0
-            else:
-                x[i, :] = state[int(j), :, 0]
-
         max_cloud_height = 3000.0
         mask[4, :] = np.tan(zen) * max_cloud_height / pixel_size
 
         # AOD 550
-        mask[5, :] = x[:, aod_bands].sum(axis=1)
-        aerosol_threshold = 0.4
+        mask[5, :] = atm[:, aod_bands].sum(axis=1)
+
+        mask[6, :] = atm[:, h2o_band].T
 
-        mask[6, :] = x[:, h2o_band].T
+        # Remove water and spacecraft flagsg if cloud flag is on (mostly cosmetic)
+        mask[2:4, np.logical_or(mask[0,:] == 1, mask[1,:] ==1)] = 0
 
-        mask[7, :] = np.array((mask[0, :] + mask[2, :] +
-                               (mask[3, :] > aerosol_threshold)) > 0, dtype=int)
         mask[:, bad] = -9999.0
         return_mask[line - start_line,...] = mask.copy()
 
@@ -128,24 +121,24 @@ def main():
     parser = argparse.ArgumentParser(description="Remove glint")
     parser.add_argument('rdnfile', type=str, metavar='RADIANCE')
     parser.add_argument('locfile', type=str, metavar='LOCATIONS')
-    parser.add_argument('lblfile', type=str, metavar='SUBSET_LABELS')
-    parser.add_argument('statefile', type=str, metavar='STATE_SUBSET')
+    parser.add_argument('atmfile', type=str, metavar='SUBSET_LABELS')
     parser.add_argument('irrfile', type=str, metavar='SOLAR_IRRADIANCE')
     parser.add_argument('outfile', type=str, metavar='OUTPUT_MASKS')
     parser.add_argument('--wavelengths', type=str, default=None)
     parser.add_argument('--n_cores', type=int, default=-1)
+    parser.add_argument('--aerosol_threshold', type=float, default=0.5)
     args = parser.parse_args()
 
     rdn_hdr = envi.read_envi_header(envi_header(args.rdnfile))
-    state_hdr = envi.read_envi_header(envi_header(args.statefile))
     rdn_shp = envi.open(envi_header(args.rdnfile)).open_memmap(interleave='bil').shape
-    lbl_shp = envi.open(envi_header(args.lblfile)).open_memmap(interleave='bil').shape
+    atm_hdr = envi.read_envi_header(envi_header(args.atmfile))
+    atm_shp = envi.open(envi_header(args.atmfile)).open_memmap(interleave='bil').shape
     loc_shp = envi.open(envi_header(args.locfile)).open_memmap(interleave='bil').shape
 
     # Check file size consistency
     if loc_shp[0] != rdn_shp[0] or loc_shp[2] != rdn_shp[2]:
         raise ValueError('LOC and input file dimensions do not match.')
-    if lbl_shp[0] != rdn_shp[0] or lbl_shp[2] != rdn_shp[2]:
+    if atm_shp[0] != rdn_shp[0] or atm_shp[2] != rdn_shp[2]:
         raise ValueError('Label and input file dimensions do not match.')
     if loc_shp[1] != 3:
         raise ValueError('LOC file should have three bands.')
@@ -165,7 +158,7 @@ def main():
 
     # Find H2O and AOD elements in state vector
     aod_bands, h2o_band = [], []
-    for i, name in enumerate(state_hdr['band names']):
+    for i, name in enumerate(atm_hdr['band names']):
         if 'H2O' in name:
             h2o_band.append(i)
         elif 'AER' in name or 'AOT' in name or 'AOD' in name:
@@ -222,26 +215,29 @@ def main():
 
     linebreaks = np.linspace(0, rdn_shp[0], num=args.n_cores*3).astype(int)
 
-    state = envi.open(envi_header(args.statefile)).open_memmap(interleave='bil').copy()
-    stateid = ray.put(state)
     irrid = ray.put(irr_resamp)
-    jobs = [build_line_masks.remote(linebreaks[_l], linebreaks[_l+1], args.rdnfile, args.locfile, args.lblfile, stateid, dt, h2o_band, aod_bands, pixel_size, args.outfile, wl, irrid) for _l in range(len(linebreaks)-1)]
+    jobs = [build_line_masks.remote(linebreaks[_l], linebreaks[_l+1], args.rdnfile, args.locfile, args.atmfile, dt, h2o_band, aod_bands, pixel_size, args.outfile, wl, irrid) for _l in range(len(linebreaks)-1)]
     rreturn = [ray.get(jid) for jid in jobs]
     ray.shutdown()
 
     mask = np.zeros((rdn_shp[0], maskbands, rdn_shp[2]))
     for lm, start_line, stop_line in rreturn:
         mask[start_line:stop_line,...] = lm
 
-    bad = np.squeeze(mask[:, 0, :]) < -9990
+    bad = np.squeeze(mask[:, 0, :]) <= -9990
     good = np.squeeze(mask[:, 0, :]) > -9990
 
-    cloudinv = np.logical_not(np.squeeze(mask[:, 0, :]))
+
+    # Create buffer around clouds (main and cirrus)
+    cloudinv = np.logical_not(np.squeeze(np.logical_or(mask[:, 0, :], mask[:,1,:])))
     cloudinv[bad] = 1
     cloud_distance = distance_transform_edt(cloudinv)
     invalid = (np.squeeze(mask[:, 4, :]) >= cloud_distance)
     mask[:, 4, :] = invalid.copy()
 
+    # Combine Cloud, Cirrus, Water, Spacecraft, and Buffer masks
+    mask[:, 7, :] = np.logical_or(np.sum(mask[:,0:5,:], axis=1) > 0, mask[:,5,:] > args.aerosol_threshold)
+
     hdr = rdn_hdr.copy()
     hdr['bands'] = str(maskbands)
     hdr['band names'] = ['Cloud flag', 'Cirrus flag', 'Water flag',

diff --git a/output_conversion.py b/output_conversion.py
@@ -22,6 +22,7 @@ def main():
     parser.add_argument('rfl_file', type=str, help="EMIT L2A reflectance ENVI file")
     parser.add_argument('rfl_unc_file', type=str, help="EMIT L2A reflectance uncertainty ENVI file")
     parser.add_argument('mask_file', type=str, help="EMIT L2A water/cloud mask ENVI file")
+    parser.add_argument('band_mask_file', type=str, help="EMIT L1B band mask ENVI file")
     parser.add_argument('loc_file', type=str, help="EMIT L1B location data ENVI file")
     parser.add_argument('glt_file', type=str, help="EMIT L1B glt ENVI file")
     parser.add_argument('version', type=str, help="3 digit (with leading V) version number")
@@ -38,6 +39,7 @@ def main():
     rfl_ds = envi.open(envi_header(args.rfl_file))
     rfl_unc_ds = envi.open(envi_header(args.rfl_unc_file))
     mask_ds = envi.open(envi_header(args.mask_file))
+    bandmask_ds = envi.open(envi_header(args.band_mask_file))
 
     # Start with Reflectance File
 
@@ -140,6 +142,7 @@ def main():
 
     logging.debug('Creating dimensions')
     makeDims(nc_ds, args.mask_file, args.glt_file)
+    nc_ds.createDimension('packed_wavelength_bands', int(bandmask_ds.metadata['bands']))
 
     logging.debug('Creating and writing mask metadata')
     add_variable(nc_ds, "sensor_band_parameters/mask_bands", str, "Mask Band Names", None,
@@ -154,6 +157,8 @@ def main():
     logging.debug('Write mask data')
     add_variable(nc_ds, 'mask', "f4", "Masks", "unitless", mask_ds.open_memmap(interleave='bip')[...].copy(),
                  {"dimensions":("downtrack", "crosstrack", "bands"), "zlib": True, "complevel": 9})
+    add_variable(nc_ds, 'band_mask', "u8", "Per-Wavelength Mask", "unitless", bandmask_ds.open_memmap(interleave='bip')[...].copy(),
+                 {"dimensions":("downtrack", "crosstrack", "packed_wavelength_bands"), "zlib": True, "complevel": 9})
     nc_ds.sync()
     nc_ds.close()
     del nc_ds

diff --git a/surface/surface_20220908.json b/surface/surface_20220908.json
@@ -0,0 +1,76 @@
+{
+  "output_model_file": "./surface.mat",
+  "wavelength_file":   "./wavelengths.txt",
+  "normalize":"Euclidean",
+  "reference_windows":[[400,1250],[1450,1750],[2050,2450]],
+  "sources":
+    [
+      {
+        "input_spectrum_files":
+          [
+            "filtered_other"
+          ],
+        "n_components": 5,
+        "windows": [
+          {"interval":[300,740], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[740,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-8, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2470,2500], "regularizer":1e-6, "correlation":"EM", "isolated": 1, "name": "noise" }
+        ]
+      },
+      {
+        "input_spectrum_files":
+          [
+            "filtered_veg"
+          ],
+        "n_components": 1,
+        "windows": [
+          {"interval":[300,350], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[350,480], "regularizer":1e-6, "correlation":"EM", "name": "aerosol"},
+          {"interval":[480,740], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[740,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-8, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2470,2500], "regularizer":1e-6, "correlation":"EM","isolated": 1, "name": "noise" }
+        ]
+      },
+      {
+        "input_spectrum_files":
+          [
+            "filtered_ocean"
+          ],
+        "n_components": 1,
+        "windows": [
+          {"interval":[300,880], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[880,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-6, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2470,2500], "regularizer":1e-6, "correlation":"EM","isolated": 1, "name": "noise" }
+        ]
+      },
+      {
+        "input_spectrum_files":
+          [
+            "surface_Liquids"
+          ],
+        "n_components": 2,
+        "windows": [
+          {"interval":[300,880], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[880,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-6, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2470,2500], "regularizer":1e-6, "correlation":"EM","isolated": 1, "name": "noise" }
+        ]
+      }
+
+  ]
+}
diff --git a/surface/surface_20221001.json b/surface/surface_20221001.json
@@ -0,0 +1,73 @@
+{
+  "output_model_file": "./surface.mat",
+  "wavelength_file":   "./wavelengths.txt",
+  "normalize":"Euclidean",
+  "reference_windows":[[400,1250],[1450,1750],[2050,2450]],
+  "sources":
+    [
+      {
+        "input_spectrum_files":
+          [
+            "filtered_other"
+          ],
+        "n_components": 5,
+        "windows": [
+          {"interval":[300,740], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[740,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-8, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2300,2500], "regularizer":1e-3, "correlation":"EM", "isolated": 1, "name": "noise" }
+        ]
+      },
+      {
+        "input_spectrum_files":
+          [
+            "filtered_veg"
+          ],
+        "n_components": 1,
+        "windows": [
+          {"interval":[300,350], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[350,480], "regularizer":1e-6, "correlation":"EM", "name": "aerosol"},
+          {"interval":[480,740], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[740,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-8, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2300,2500], "regularizer":1e-3, "correlation":"EM","isolated": 1, "name": "noise" }
+        ]
+      },
+      {
+        "input_spectrum_files":
+          [
+            "filtered_ocean"
+          ],
+        "n_components": 1,
+        "windows": [
+          {"interval":[300,740], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[740,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-6, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,2500], "regularizer":10, "correlation": "decorrelated" }
+        ]
+      },
+      {
+        "input_spectrum_files":
+          [
+            "surface_Liquids"
+          ],
+        "n_components": 2,
+        "windows": [
+          {"interval":[300,740], "regularizer":10, "correlation":"decorrelated"},
+          {"interval":[740,1250], "regularizer":1e-6, "correlation":"EM", "name": "shallow-water"},
+          {"interval":[1250,1325], "regularizer":1e-6, "correlation":"EM", "name": "osf"},
+          {"interval":[1325,1960], "regularizer":10, "correlation": "decorrelated" },
+          {"interval":[1960,2070], "regularizer":1e-6, "correlation":"EM","name": "co2" },
+          {"interval":[2070,2470], "regularizer":10, "correlation":"decorrelated" },
+          {"interval":[2300,2500], "regularizer":1e-3, "correlation":"EM","isolated": 1, "name": "noise" }
+        ]
+      }
+
+  ]
+}