still working on the scene refactoring

mvcm/constants.py

@@ -34,6 +34,14 @@ class SceneConstants:
     POLAR_LATITUDES = 60.0
     POLAR_COAST_LATITUDES = 64.0
 
+    ARCTIC_LAT = 67.5
+    GREENLAND_BB_W_LON = -40.0
+    GREENLAND_BB_E_LON = -30.0
+    SIBERIA_BB_W_LON = -12.5
+    CANADA_BB_W_LON = -168.5
+
+    VIS_RATIO = 0.9
+
     SCENE_LIST = (
         "Ocean_Day",
         "Ocean_Night",

mvcm/scene.py

@@ -139,8 +139,8 @@ def test_scene():
 """
 
 
-def compute_ref_angle(vza, sza, raz):
-    """Compute ref angle."""
+def compute_glint_angle(vza, sza, raz):
+    """Compute glint angle."""
     cos_refang = np.sin(vza * Constants.DTR) * np.sin(sza * Constants.DTR) * np.cos(
         raz * Constants.DTR
     ) + np.cos(vza * Constants.DTR) * np.cos(sza * Constants.DTR)
@@ -158,7 +158,7 @@ class Scene:
     ref_ang: np.ndarray = field(
         init=False,
         default=Factory(
-            lambda self: compute_ref_angle(
+            lambda self: compute_glint_angle(
                 self.data.sensor_zenith.values,
                 self.data.solar_zenith.values,
                 self.data.relative_azimuth.values,
@@ -354,10 +354,27 @@ def find_scene(
         & (data.land_water_mask != land_water_flags["ephemeral"]),
         1,
         scene_flag["land"],
+    )
+
+    scene_flag["land"] = xr.where(
+        (data.land_water_mask == SceneConstants.COASTLINE)
+        | (data.land_water_mask == SceneConstants.SHALLOW_INLAND),
+        1,
+        scene_flag["land"],
+    )
+
+    scene_flag["land"] = xr.where(
+        (data.land_water_mask == SceneConstants.NO_DATA)
+        & (data.M07 != Constants.BAD_DATA)
+        & (data.M05 != Constants.BAD_DATA)
+        & (data.M07 / data.M05 > SceneConstants.VIS_RATIO),
+        1,
+        0,
     ).values
 
     # idx = np.nonzero(((lsf == 1) | (lsf == 4)) & (eco == 14) & (lat < 64.0))
     # scene_flag["coast"][idx] = 1
+    # Coast for everything not in polar regiones
     scene_flag["coast"] = xr.where(
         (
             (data.land_water_mask == land_water_flags["land"])
@@ -373,37 +390,69 @@ def find_scene(
     #     (eco == 14) & ((lsf == 1) | (lsf == 4)) & ((lat >= 67.5) & (lon < -40.0) & (lon > -168.6))
     # )
     # scene_flag["coast"][idx] = 1
+
+    # idx = np.nonzero((eco == 14) & ((lsf == 1) | (lsf == 4)) & ((lat >= 67.5) & (lon > -12.5)))
+    # scene_flag["coast"][idx] = 1
+    # idx = np.nonzero(
+    #     (eco == 14)
+    #     & ((lsf == 1) | (lsf == 4))
+    #     & ((lat >= 64.0) & (lat < 67.5) & (lon < -40.0) & (lon > -168.5))
+    # )
+    # scene_flag["coast"][idx] = 1
+    # idx = np.nonzero(
+    #     (eco == 14) & ((lsf == 1) | (lsf == 4)) & ((lat >= 64.0) & (lat < 67.5) & (lon > -30.0))
+    # )
+    # scene_flag["coast"][idx] = 1
+
+    # The colors in the comments below refer to the highlighted areas as shown in the map provided
+    # in the documentation (figures/coast_areas_scene_selection_mvcm). The figure is used to
+    # remember the logic of this block that otherwise is rather confusing
     scene_flag["coast"] = xr.where(
         (
-            (data.ecosystem == ecosystem_flags["water"])
-            & (data.land_water_mask == land_water_flags["land"])
-            | (data.land_water_mask == land_water_flags["ephemeral"])
+            (data.ecosystem == SceneConstants.ECO_FLAG_WATER)
+            & (data.land_water_mask == SceneConstants.LAND)
+            | (data.land_water_mask == SceneConstants.EPHEMERAL)
+        )
+        & (  # GREEN
+            (data.latitude >= SceneConstants.ARCTIC_LAT)
+            & (
+                (data.longitude < SceneConstants.GREENLAND_BB_W_LON)
+                & (data.longitude > SceneConstants.CANADA_BB_W_LON)
+            )
+            | (data.longitude > SceneConstants.SIBERIA_BB_W_LON)
         )
-        & ((data.latitude >= 67.5) & (data.longitude < -40.0) & (data.longitude > -168.6)),
+        | (  # YELLOW
+            (data.latitude >= SceneConstants.POLAR_COAST_LATITUDES)
+            & (data.latitude < SceneConstants.ARCTIC_LAT)
+            & (  # RED
+                (data.longitude < SceneConstants.GREENLAND_BB_W_LON)
+                & (data.longitude > SceneConstants.CANADA_BB_W_LON)
+                | (data.longitude > SceneConstants.GREENLAND_BB_E_LON)
+            )
+        ),
+        1,
+        scene_flag["coast"],
+    )
+
+    scene_flag["coast"] = xr.where(
+        (data.land_water_mask == SceneConstants.COASTLINE)
+        | ((data.land_water_mask == SceneConstants.SHALLOW_INLAND) & (scene_flag["day"] == 1)),
         1,
         scene_flag["coast"],
     ).values
 
-    idx = np.nonzero((eco == 14) & ((lsf == 1) | (lsf == 4)) & ((lat >= 67.5) & (lon > -12.5)))
-    scene_flag["coast"][idx] = 1
+    # ARCTIC_LAT = 67.5
+    # GREENLAND_BB_W_LON = -40.0
+    # GREENLAND_BB_E_LON = -30.0
+    # SIBERIA_BB_W_LON = -12.5
+    # CANADA_BB_W_LON = -168.5
 
-    idx = np.nonzero(
-        (eco == 14)
-        & ((lsf == 1) | (lsf == 4))
-        & ((lat >= 64.0) & (lat < 67.5) & (lon < -40.0) & (lon > -168.5))
-    )
-    scene_flag["coast"][idx] = 1
-    idx = np.nonzero(
-        (eco == 14) & ((lsf == 1) | (lsf == 4)) & ((lat >= 64.0) & (lat < 67.5) & (lon > -30.0))
-    )
-    scene_flag["coast"][idx] = 1
-
-    idx = np.nonzero(lsf == 2)
-    scene_flag["coast"][idx] = 1
-    scene_flag["land"][idx] = 1
+    # idx = np.nonzero(lsf == 2)
+    # scene_flag["coast"][idx] = 1
+    # scene_flag["land"][idx] = 1
 
     idx = np.nonzero(lsf == 3)
-    scene_flag["land"][idx] = 1
+    # scene_flag["land"][idx] = 1
     scene_flag["sh_lake"][idx] = 1
 
     idx = np.nonzero((lsf == 0) | ((lsf >= 5) & (lsf <= 7)))
@@ -413,18 +462,18 @@ def find_scene(
     scene_flag["sh_ocean"][lsf == 0] = 1
 
     # Need shallow lakes to be processed as "coast" for day, but not night
-    idx = np.nonzero((lsf == 3) & (day == 1))
-    scene_flag["coast"][idx] = 1
+    # idx = np.nonzero((lsf == 3) & (day == 1))
+    # scene_flag["coast"][idx] = 1
 
     # if land/sea flag is missing, then calculate visible ratio
     # to determine if land or water.
-    idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086 / b065 > 0.9))
-    scene_flag["land"][idx] = 1
+    # idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086 / b065 > 0.9))
+    # scene_flag["land"][idx] = 1
 
     idx = np.nonzero(
         (lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086 / b065 <= 0.9)
     )
-    scene_flag["land"][idx] = 0
+    # scene_flag["land"][idx] = 0 # I implemented this, but not the water part
     scene_flag["water"][idx] = 1
 
     # Check surface elevation