Merge pull request #56 from Open-ET/et-fraction-grass-source-param

Add support for an ET fraction grass source parameter

Author: cgmorton

openet/ssebop/image.py

@@ -51,6 +51,7 @@ def __init__(
             dt_min=5,
             dt_max=25,
             et_fraction_type='alfalfa',
+            et_fraction_grass_source=None,
             reflectance_type='SR',
             **kwargs,
     ):
@@ -85,7 +86,7 @@ def __init__(
                        'DAYMET_MEDIAN_V2', 'CIMIS_MEDIAN_V1',
                        collection ID, or float}, optional
             Maximum air temperature source.  The default is
-            'projects/usgs-ssebop/tmax/daymet_v3_median_1980_2018'.
+            'projects/usgs-ssebop/tmax/daymet_v4_mean_1981_2010'.
         elr_flag : bool, str, optional
             If True, apply Elevation Lapse Rate (ELR) adjustment
             (the default is False).
@@ -94,10 +95,16 @@ def __init__(
         dt_max : float, optional
             Maximum allowable dT [K] (the default is 25).
         et_fraction_type : {'alfalfa', 'grass'}, optional
-            ET fraction  (the default is 'alfalfa').
+            ET fraction reference type (the default is 'alfalfa').
+            If set to "grass" the et_fraction_grass_source must also be set.
+        et_fraction_grass_source : {'NASA/NLDAS/FORA0125_H002',
+                                    'ECMWF/ERA5_LAND/HOURLY'}, float, optional
+            Reference ET source for alfalfa to grass reference adjustment.
+            Parameter must be set if et_fraction_type is 'grass'.
+            The default is currently the NLDAS hourly collection,
+            but having a default will likely be removed in a future version.
         reflectance_type : {'SR', 'TOA'}, optional
-            Used to set the Tcorr NDVI thresholds
-            (the default is 'SR').
+            Used to set the Tcorr NDVI thresholds (the default is 'SR').
         kwargs : dict, optional
             dt_resample : {'nearest', 'bilinear'}
             tcorr_resample : {'nearest', 'bilinear'}
@@ -194,14 +201,28 @@ def __init__(
                 raise ValueError(f'elr_flag "{self._elr_flag}" could not be interpreted as bool')
 
         # ET fraction type
-        # CGM - Should et_fraction_type be set as a kwarg instead?
         if et_fraction_type.lower() not in ['alfalfa', 'grass']:
             raise ValueError('et_fraction_type must "alfalfa" or "grass"')
         self.et_fraction_type = et_fraction_type.lower()
-        # if 'et_fraction_type' in kwargs.keys():
-        #     self.et_fraction_type = kwargs['et_fraction_type'].lower()
-        # else:
-        #     self.et_fraction_type = 'alfalfa'
+
+        # ET fraction alfalfa to grass reference adjustment
+        # The NLDAS hourly collection will be used if a source value is not set
+        if self.et_fraction_type.lower() == 'grass' and not et_fraction_grass_source:
+            warnings.warn(
+                'NLDAS is being set as the default ET fraction grass adjustment source.  '
+                'In a future version the parameter will need to be set explicitly as: '
+                'et_fraction_grass_source="NASA/NLDAS/FORA0125_H002".',
+                FutureWarning
+            )
+            et_fraction_grass_source = 'NASA/NLDAS/FORA0125_H002'
+        self.et_fraction_grass_source = et_fraction_grass_source
+        # if self.et_fraction_type.lower() == 'grass' and not et_fraction_grass_source:
+        #     raise ValueError(
+        #         'et_fraction_grass_source parameter must be set if et_fraction_type==\'grass\''
+        #     )
+        # # Should the supported source values be checked here instead of in model.py?
+        # if et_fraction_grass_source not in et_fraction_grass_sources:
+        #     raise ValueError('unsupported et_fraction_grass_source')
 
         self.reflectance_type = reflectance_type
         if reflectance_type not in ['SR', 'TOA']:
@@ -333,50 +354,16 @@ def et_fraction(self):
 
         et_fraction = model.et_fraction(lst=self.lst, tmax=tmax, tcorr=self.tcorr, dt=dt)
 
-        # TODO: Add support for setting the conversion source dataset
-        # TODO: Interpolate "instantaneous" ETo and ETr?
-        # TODO: Move openet.refetgee import to top?
-        # TODO: Check if etr/eto is right (I think it is)
-        if self.et_fraction_type.lower() == 'grass':
-            import openet.refetgee
-            nldas_coll = (
-                ee.ImageCollection('NASA/NLDAS/FORA0125_H002')
-                .select(['temperature', 'specific_humidity', 'shortwave_radiation', 'wind_u', 'wind_v'])
-            )
-
-            # Interpolating hourly NLDAS to the Landsat scene time
-            # CGM - The 2 hour window is useful in case an image is missing
-            #   I think EEMETRIC is using a 4 hour window
-            # CGM - Need to check if the NLDAS images are instantaneous
-            #   or some sort of average of the previous or next hour
-            time_start = ee.Number(self._time_start)
-            prev_img = ee.Image(
-                nldas_coll
-                .filterDate(time_start.subtract(2 * 60 * 60 * 1000), time_start)
-                .limit(1, 'system:time_start', False)
-                .first()
-            )
-            next_img = ee.Image(
-                nldas_coll.filterDate(time_start, time_start.add(2 * 60 * 60 * 1000)).first()
-            )
-            prev_time = ee.Number(prev_img.get('system:time_start'))
-            next_time = ee.Number(next_img.get('system:time_start'))
-            time_ratio = time_start.subtract(prev_time).divide(next_time.subtract(prev_time))
-            nldas_img = (
-                next_img.subtract(prev_img).multiply(time_ratio).add(prev_img)
-                .set({'system:time_start': self._time_start})
-            )
-
-            # # DEADBEEF - Select NLDAS image before the Landsat scene time
-            # nldas_img = ee.Image(nldas_coll
-            #     .filterDate(self._date.advance(-1, 'hour'), self._date)
-            #     .first())
-
-            et_fraction = (
-                et_fraction
-                .multiply(openet.refetgee.Hourly.nldas(nldas_img).etr)
-                .divide(openet.refetgee.Hourly.nldas(nldas_img).eto)
-            )
+        # Convert the ET fraction to a grass reference fraction
+        if self.et_fraction_type.lower() == 'grass' and self.et_fraction_grass_source:
+            if utils.is_number(self.et_fraction_grass_source):
+                et_fraction = et_fraction.multiply(self.et_fraction_grass_source)
+            else:
+                et_fraction = model.etf_grass_type_adjust(
+                    etf=et_fraction,
+                    src_coll_id=self.et_fraction_grass_source,
+                    time_start=self._time_start
+                )
 
         return et_fraction.set(self._properties)\
             .set({

openet/ssebop/model.py

@@ -2,6 +2,8 @@
 
 import ee
 
+import openet.refetgee
+
 
 def et_fraction(lst, tmax, tcorr, dt):
     """SSEBop fraction of reference ET (ETf)
@@ -31,17 +33,12 @@ def et_fraction(lst, tmax, tcorr, dt):
 
     """
 
-    et_fraction = lst.expression(
+    etf = lst.expression(
         '(lst * (-1) + tmax * tcorr + dt) / dt',
         {'tmax': tmax, 'dt': dt, 'lst': lst, 'tcorr': tcorr}
     )
 
-    return (
-        et_fraction
-        .updateMask(et_fraction.lte(2.0))
-        .clamp(0, 1.0)
-        .rename(['et_fraction'])
-    )
+    return etf.updateMask(etf.lte(2.0)).clamp(0, 1.0).rename(['et_fraction'])
 
 
 def dt(tmax, tmin, elev, doy, lat=None, rs=None, ea=None):
@@ -147,9 +144,7 @@ def dt(tmax, tmin, elev, doy, lat=None, rs=None, ea=None):
     den = tmax.add(tmin).multiply(0.5).pow(-1).multiply(pair).multiply(3.486 / 1.01)
 
     # Temperature difference [K] (Senay2018 A.5)
-    dt = rn.divide(den).multiply(110.0 / ((1.013 / 1000) * 86400))
-
-    return dt
+    return rn.divide(den).multiply(110.0 / ((1.013 / 1000) * 86400))
 
 
 def lapse_adjust(temperature, elev, lapse_threshold=1500):
@@ -191,7 +186,7 @@ def elr_adjust(temperature, elevation, radius=80):
 
     Returns
     -------
-    ee.Image of adjusted temperature
+    ee.Image
 
     Notes
     -----
@@ -233,3 +228,79 @@ def elr_adjust(temperature, elevation, radius=80):
     tmax_img = tmax_img.where(elr_mask, elr_adjust)
 
     return tmax_img
+
+
+# TODO: Decide if using the instantaneous is the right/best approach
+#   We could use the closest hour in time, an average of a few hours
+#   or just switch to using the raw daily or bias corrected assets
+def etf_grass_type_adjust(etf, src_coll_id, time_start):
+    """"Convert ET fraction from an alfalfa reference to grass reference
+
+    Parameters
+    ----------
+    etf : ee.Image
+        ET fraction (alfalfa reference).
+    src_coll_id : str
+        Hourly meteorology collection ID for computing reference ET.
+    time_start : int, ee.Number
+        Image system time start [millis].
+
+    Returns
+    -------
+    ee.Image
+
+    """
+    hourly_et_reference_sources = [
+        'NASA/NLDAS/FORA0125_H002',
+        'ECMWF/ERA5_LAND/HOURLY',
+    ]
+    if src_coll_id not in hourly_et_reference_sources:
+        raise ValueError(f'unsupported hourly ET reference source: {src_coll_id}')
+    elif not src_coll_id:
+        raise ValueError('hourly ET reference source not')
+    else:
+        src_coll = ee.ImageCollection(src_coll_id)
+
+    # Interpolating hourly NLDAS to the Landsat scene time
+    # CGM - The 2 hour window is useful in case an image is missing
+    #   I think EEMETRIC is using a 4 hour window
+    # CGM - Need to check if the NLDAS images are instantaneous
+    #   or some sort of average of the previous or next hour
+    time_start = ee.Number(time_start)
+    prev_img = ee.Image(
+        src_coll
+        .filterDate(time_start.subtract(2 * 60 * 60 * 1000), time_start)
+        .limit(1, 'system:time_start', False)
+        .first()
+    )
+    next_img = ee.Image(
+        src_coll.filterDate(time_start, time_start.add(2 * 60 * 60 * 1000)).first()
+    )
+    prev_time = ee.Number(prev_img.get('system:time_start'))
+    next_time = ee.Number(next_img.get('system:time_start'))
+    time_ratio = time_start.subtract(prev_time).divide(next_time.subtract(prev_time))
+    interp_img = (
+        next_img.subtract(prev_img).multiply(time_ratio).add(prev_img)
+        .set({'system:time_start': time_start})
+    )
+
+    # # DEADBEEF - Select the NLDAS image before the Landsat scene time
+    # interp_img = ee.Image(
+    #     hourly_coll.filterDate(self._date.advance(-1, 'hour'), self._date).first()
+    # )
+
+    if src_coll_id.upper() == 'NASA/NLDAS/FORA0125_H002':
+        etf_grass = (
+            etf
+            .multiply(openet.refetgee.Hourly.nldas(interp_img).etr)
+            .divide(openet.refetgee.Hourly.nldas(interp_img).eto)
+        )
+    elif src_coll_id.upper() == 'ECMWF/ERA5_LAND/HOURLY':
+        etf_grass = (
+            etf
+            .multiply(openet.refetgee.Hourly.era5_land(interp_img).etr)
+            .divide(openet.refetgee.Hourly.era5_land(interp_img).eto)
+        )
+    # else:
+
+    return etf_grass

openet/ssebop/tests/test_b_model.py

@@ -4,6 +4,16 @@
 import openet.ssebop.model as model
 import openet.ssebop.utils as utils
 
+COLL_ID = 'LANDSAT/LC08/C02/T1_L2'
+SCENE_ID = 'LC08_042035_20150713'
+# SCENE_DT = datetime.datetime.strptime(SCENE_ID[-8:], '%Y%m%d')
+# SCENE_DATE = SCENE_DT.strftime('%Y-%m-%d')
+# SCENE_DOY = int(SCENE_DT.strftime('%j'))
+# # SCENE_TIME = utils.millis(SCENE_DT)
+SCENE_TIME = 1436812419150
+SCENE_POINT = (-119.5, 36.0)
+TEST_POINT = (-119.44252382373145, 36.04047742246546)
+
 
 @pytest.mark.parametrize(
     # Note: These are made up values
@@ -187,3 +197,48 @@ def test_Model_elr_adjust(xy, adjusted):
         assert output < original
     else:
         assert output == original
+
+
+def test_Image_et_reference_source_parameters():
+    """Check that the function parameter names and order don't change"""
+    etf_img = (
+        ee.Image(f'{COLL_ID}/{SCENE_ID}').select([0]).multiply(0).add(1.0)
+        .rename(['et_fraction']).set('system:time_start', SCENE_TIME)
+    )
+    output = model.etf_grass_type_adjust(
+        etf=etf_img, src_coll_id='NASA/NLDAS/FORA0125_H002', time_start=SCENE_TIME
+    )
+    assert utils.point_image_value(output, SCENE_POINT, scale=100)['et_fraction'] > 1
+
+    output = model.etf_grass_type_adjust(etf_img, 'NASA/NLDAS/FORA0125_H002', SCENE_TIME)
+    assert utils.point_image_value(output, SCENE_POINT, scale=100)['et_fraction'] > 1
+
+
+@pytest.mark.parametrize(
+    'src_coll_id, expected',
+    [
+        ['NASA/NLDAS/FORA0125_H002', 1.23],
+        ['ECMWF/ERA5_LAND/HOURLY', 1.15],
+    ]
+)
+def test_Model_etf_grass_type_adjust(src_coll_id, expected, tol=0.01):
+    """Check alfalfa to grass reference adjustment factor"""
+    etf_img = (
+        ee.Image(f'{COLL_ID}/{SCENE_ID}').select([0]).multiply(0).add(1.0)
+        .rename(['et_fraction']).set('system:time_start', SCENE_TIME)
+    )
+    output = model.etf_grass_type_adjust(
+        etf=etf_img, src_coll_id=src_coll_id, time_start=SCENE_TIME
+    )
+    output = utils.point_image_value(output, SCENE_POINT, scale=100)
+    assert abs(output['et_fraction'] - expected) <= tol
+
+
+def test_Model_etf_grass_type_adjust_src_coll_id_exception():
+    """Function should raise an exception for unsupported src_coll_id values"""
+    with pytest.raises(ValueError):
+        utils.getinfo(model.etf_grass_type_adjust(
+            etf=ee.Image.constant(1), src_coll_id='DEADBEEF', time_start=SCENE_TIME
+        ))
+
+