Merge pull request #13 from lopezvoliver/eelandsat

ET estimation workflow with optional vegetation and positive LE masking

geeet/eepredefined/__init__.py

@@ -1,4 +1,6 @@
 import geeet.eepredefined.landsat
 import geeet.eepredefined.join
+import geeet.eepredefined.pixel_area
+import geeet.eepredefined.masks
 from geeet.eepredefined.metprep import MeteoBands, MeteoPrep
 

geeet/eepredefined/landsat.py

@@ -252,6 +252,7 @@ def collection(
     Returns: ee.ImageCollection
     """
     import ee
+    from .join import landsat_ecmwf
     from .parsers import feature_collection
 
     region = feature_collection(region)
@@ -337,7 +338,7 @@ def collection(
             .select(*meteo_bands)
             .map(meteo_prep))
 
-        collection = geeet.eepredefined.join.landsat_ecmwf(
+        collection = landsat_ecmwf(
             collection, meteo_collection)
 
         # Set ERA5 measuring heights (zU, zT)

geeet/eepredefined/masks.py

@@ -0,0 +1,67 @@
+"""
+Custom masks (based on NDVI, positive LE, fractional vegetation cover)
+"""
+import ee
+from typing import Callable
+
+def apply_static_mask(mask:str, bandNames:list)->Callable:
+    def apply_mask(img:ee.Image)->ee.Image:
+        """Returns img with the given mask (str) used to 
+        update the mask on selected bands (bandNames)
+        Requires bands: [mask]
+        """
+        imgBands = (img.select(bandNames)
+        .updateMask(img.select(mask))
+        )
+        return img.addBands(imgBands, overwrite=True)
+
+    return apply_mask
+
+
+def Fndvi_mask(NDVI_BARE_GROUND=0.2)->Callable:
+    def ndvi_mask(img:ee.Image)->ee.Image:
+        """Returns a mask based on a threshold for NDVI
+        Requires bands: "NDVI"
+        Adds band: "vegetation_mask"
+        """
+        return img.addBands(
+           img.select("NDVI").gt(NDVI_BARE_GROUND) 
+           .rename("vegetation_mask")
+        )
+    return ndvi_mask
+
+
+def Ffvc(fvc:ee.ImageCollection)->Callable:
+    """Returns a mappable function to add fractional vegetation cover (fvc)
+    
+    Given an image collection of a single-band image (fvc), 
+    this function adds fvc to an image for the given year, 
+    **only** within the **unobserved** regions. 
+
+    Requires band: "cloud_cover"
+    
+    Adds band: "fvc"
+    """
+    fvc = ee.ImageCollection(fvc)
+    def f(img:ee.Image)->ee.Image:
+        year=ee.Number(img.date().get("year"))
+        return img.addBands(
+            fvc.filter(ee.Filter.calendarRange(year, year, "year"))
+            .first()
+            .updateMask(
+                img.select("cloud_cover")  
+                .unmask(1)
+            )
+            .rename("fvc") 
+        )
+    return f
+
+
+def positive_LE_mask(img:ee.Image)->ee.Image:
+    """
+    Mask to keep only positive LE values. 
+    """
+    return img.addBands(
+        img.select("LE").gt(0)
+        .rename("positive_le_mask")
+    )

geeet/eepredefined/pixel_area.py

@@ -0,0 +1,70 @@
+"""
+Custom pixel area mappable functions. 
+"""
+import ee
+
+def feature_area(img: ee.Image)->ee.Image:
+    """Returns feature_area (pixelArea) as a band
+    Reduce this band using ee.Reducer.sum() to get 
+    the total area in m².
+
+    Requires bands: None
+    Adds bands: "feature_area"
+    """
+    return img.addBands(
+        ee.Image.pixelArea().rename("feature_area")
+    )
+
+
+def unobserved_area(img: ee.Image)->ee.Image:
+    """Returns unobserved_area (pixelArea) as a band
+    
+    For Landsat images, unobserved area is due to 
+    cloud/cloud shadow mask and Landsat 7 slc error 
+    gaps (stripes). 
+
+    Reduce this band using ee.Reducer.sum() to get 
+    the total unobserved area in m².
+    
+    Requires band: "cloud_cover", "feature_area"
+    
+    Adds band: "unobserved"
+    """
+    feature_pixel_area = img.select("feature_area")
+    return img.addBands(
+        feature_pixel_area.updateMask(
+            img.select("cloud_cover")
+            .unmask(1)
+        )
+        .rename("unobserved_area")
+    )
+
+
+def observed_veg_area(img:ee.Image)->ee.Image:
+    """Returns observed_vegetation_area (pixelArea) as a band
+
+    Requires band: "vegetation_mask", "feature_area"
+    
+    Adds band: "observed_vegetation_area"
+    """
+    feature_pixel_area = img.select("feature_area")
+    return img.addBands(
+        feature_pixel_area.updateMask(
+            img.select("vegetation_mask")
+        )
+        .rename("observed_vegetation_area")
+    )
+
+def unobserved_veg_area(img:ee.Image)->ee.Image:
+    """Adds the unobserved vegetation area band
+    Requires band: "fvc", "feature_area"
+    Adds band: "unobserved_vegetation_area"
+    """
+    feature_pixel_area = img.select("feature_area")
+    return img.addBands(
+        feature_pixel_area.updateMask(
+            img.select("fvc")
+        )
+        .rename("unobserved_vegetation_area")
+    )
+

geeet/eepredefined/workflows.py

@@ -0,0 +1,64 @@
+"""
+Custom workflows 
+"""
+
+def masked_et(cfmaskable_bands,
+               maskable_bands,
+               positive_le = True,
+               NDVI_BARE_GROUND = None,
+               fvc = None
+               ):
+    """
+    Custom ET estimation workflow
+
+    1. Extrapolate LE (W/m²) -> ET (mm/day)
+    2. Applies cloud-mask to cfmaskable_bands
+    3. Adds the "feature_area" pixelArea band.
+    4. Adds the "unobserved_area" pixelArea band.
+
+    If positive_le is set to True:
+        5. Adds a "positive_le_mask"
+        6. Applies the "positive_le_mask" to maskable_bands
+    
+    If NDVI_BARE_GROUND is not None:
+        7. Adds "vegetation_mask"
+        8. Applies the "vegetation_mask" to maskable_bands
+        9. Adds the "observed_vegetation_area" band.
+        
+    If fvc is not None:
+        10. Adds the fractional vegetation cover (fvc) band. 
+        11. Adds the "unobserved_vegetation_area" band. 
+    """
+    from . import landsat
+    from . import pixel_area
+    from . import masks
+
+    cfmask = landsat.cfmask(cfmaskable_bands)
+    lemask = masks.apply_static_mask("positive_le_mask", maskable_bands)
+    vegmask = masks.apply_static_mask("vegetation_mask", maskable_bands)
+    add_veg_mask = masks.Fndvi_mask(NDVI_BARE_GROUND)
+    add_fvc = masks.Ffvc(fvc)
+
+    workflow = [landsat.extrapolate_LE,
+                cfmask,
+                pixel_area.feature_area,
+                pixel_area.unobserved_area
+                ]
+
+    if positive_le:
+        workflow = workflow + [
+            masks.positive_LE_mask, lemask
+        ]
+
+    if NDVI_BARE_GROUND:
+        workflow = workflow + [
+            add_veg_mask, vegmask,
+            pixel_area.observed_veg_area
+        ]    
+
+    if fvc:
+        workflow = workflow + [
+            add_fvc, pixel_area.unobserved_veg_area, 
+        ]
+
+    return workflow