Merge pull request #25 from ikramelhazdour/meteo

geeet/eepredefined/metprep.py

@@ -3,10 +3,15 @@
 functions for mapping to an Image Collection
 """
 class MeteoBands:
+    """
+    Class to rename bands from 
+    https://developers.google.com/earth-engine/datasets/catalog/ECMWF_ERA5_LAND_HOURLY#bands
+    """
     ECMWF_ERA5_HOURLY_TSEB=(
     [
     'surface_pressure', 
     'temperature_2m', 
+    'dewpoint_temperature_2m',
     'u_component_of_wind_10m',
     'v_component_of_wind_10m',
     'surface_solar_radiation_downwards_hourly',
@@ -15,6 +20,7 @@ class MeteoBands:
     [
     'surface_pressure',
     'air_temperature', 
+    'dewpoint_temperature',
     'u_component_of_wind_10m',
     'v_component_of_wind_10m',
     'surface_solar_radiation_downwards_hourly',

geeet/meteo.py

@@ -19,10 +19,17 @@ def teten(T):
     '''
     Compute Teten's formula for saturation water vapour pressure  (esat (T)) in Pa 
     with parameters set according to Buck (1981) for saturation over water. 
-    Reference: 
-    https://www.ecmwf.int/sites/default/files/elibrary/2016/17117-part-iv-physical-processes.pdf#subsection.7.2.1    
+
+    Note that dewpoint temperature (Td) is defined as the temperature at which a sample of air
+    with vapor pressure *e* must be cooled to become saturated, 
+        e = esat(Td)
+    Therefore this formula can be used to get actual water vapor pressure (ea) from Td. 
 
     Input: T (numpy array or ee.Image) Temperature in Kelvin
+
+    References:
+    https://www.ecmwf.int/sites/default/files/elibrary/2016/17117-part-iv-physical-processes.pdf#subsection.7.2.1 
+    https://doi.org/10.1016/B978-0-12-386910-4.00002-0    
     '''
     if is_img(T):
         T1 = T.subtract(T0) # in K
@@ -32,44 +39,45 @@ def teten(T):
         esat = a1*np.exp(a3*(T-T0)/(T-a4))
     return esat
 
-def specific_humidity(T,P):
+def specific_humidity(esat,P):
     '''
+    Saturation specific humidity
+
     Input: (ee.Images or np.arrays):
+        - esat: saturation water vapor pressure
         - P: surface pressure in Pascals
-        - T: temperature in Kelvin 
     Output: (ee.Image or np.array):
         - Q: specific humidity 
+
+    References:
+    https://www.ecmwf.int/sites/default/files/elibrary/2016/17117-part-iv-physical-processes.pdf#subsection.7.2.1 
     '''
-    if is_img(T):
-        esat = teten(T)
+    if is_img(esat):
         denom = P.subtract(esat.multiply(1-epsilon))
         Q = esat.multiply(epsilon).divide(denom)
     else:
-        esat = teten(T)
         Q = epsilon*esat/(P-(1-epsilon)*esat)
     return Q
 
-def relative_humidity(temperature, dewpoint_temperature, pressure, band_name='relative_humidity'):
+def relative_humidity(temperature, dewpoint_temperature, band_name='relative_humidity'):
     '''
     Input: (ee.Images or np.arrays):
         - temperature, in Kelvin
         - dewpoint_temperature, in Kelvin
-        - pressure: surface pressure in Pascals
     Output: (ee.Image or np.array):
         - RH: relative humidity(%)
     
     Equation 7.91 in:
     https://www.ecmwf.int/sites/default/files/elibrary/2016/17117-part-iv-physical-processes.pdf
     '''
     if is_img(temperature):
-        Q = specific_humidity(dewpoint_temperature,pressure)
+        ea = teten(dewpoint_temperature)
         esat = teten(temperature)
-        denom = Q.multiply(1/epsilon -1).add(1).multiply(esat)
-        RH = pressure.multiply(Q).multiply(100/epsilon).divide(denom).rename(band_name)
+        RH = ea.multiply(100).divide(esat).rename(band_name)
     else:
+        ea = teten(dewpoint_temperature)
         esat = teten(temperature)
-        Q = specific_humidity(dewpoint_temperature,pressure)
-        RH = (pressure*Q*100/epsilon)/(esat*(1+Q*((1/epsilon) - 1)))
+        RH = 100*ea/esat
     return RH
 
 def vpd(RH, Temp_K, band_name=None):
@@ -124,14 +132,15 @@ def LatHeatVap(Temp_K):
         L = 2.501 - (2.361e-3*(Temp_K-273.15)) # at 20C this is ~2.45 MJ kg-1
     return L
 
-def compute_met_params(temperature, pressure):
+def compute_met_params(temperature, dewpoint_temperature, pressure):
     """
     Calculates several temperature and/or pressure-dependent
     parameters related to heat flux in air,
     which are commonly used in ET models
 
     Inputs (ee.Image or np.arrays):
-    - temperature: air temperature at reference height (Kelvin).
+    - temperature: air temperature at reference height (Kelvin)
+    - dewpoint_temperature: Dewpoint temperature at reference height (Kelvin)
     - pressure: total air pressure (dry air + water vapor) (Pa)
     
     Outputs: (ee.Image with following bands, OR list of np.arrays:)
@@ -143,9 +152,13 @@ def compute_met_params(temperature, pressure):
     - lambda (latent heat of vaporization), in MJ kg-1
     - psicr (psicrometric constant), in Pa K-1
     - taylor (=s/(s+psicr)), in Pa K-1
+
+    References
+    - Hydrology - An Introduction (Brutsaert 2005)
+    - https://www.ecmwf.int/sites/default/files/elibrary/2016/17117-part-iv-physical-processes.pdf#section.2.7 
     """
-    q = specific_humidity(temperature, pressure) 
-    ea = teten(temperature)  # in Pa
+    ea = teten(dewpoint_temperature)  # in Pa
+    q = specific_humidity(ea, pressure)
     Lambda = LatHeatVap(temperature)  # in MJ kg-1
 
     if is_img(pressure):

geeet/tseb.py

@@ -93,7 +93,7 @@ def update_var(var, pixelsToUpdate, updated_var):
 
 def tseb_series(img=None,    # ee.Image with inputs as bands (takes precedence over numpy arrays)
                 Tr=None, NDVI=None, LAI = None,  # numpy arrays
-                P = None, Ta = None, U = None, Sdn = None, Ldn = None, Rn = None,
+                P = None, Ta = None, Td = None, U = None, Sdn = None, Ldn = None, Rn = None,
                 Alb = None,
                 doy = None, time = None, Vza = None, longitude = None, latitude=None,
                 CH = 0.3, F_g = 1.0, k_par = 0.5, k_rns=0.45,    
@@ -186,6 +186,7 @@ def tseb_series(img=None,    # ee.Image with inputs as bands (takes precedence o
         NDVI = img.select('NDVI')
         Tr = img.select('radiometric_temperature') # in K
         Ta = img.select('air_temperature')   # in K
+        Td = img.select('dewpoint_temperature')   # in K
         P = img.select('surface_pressure')   # in Pa
         U = img.select('wind_speed')         # in m/s
         Sdn = img.select('solar_radiation')  # in W/m2
@@ -210,6 +211,7 @@ def tseb_series(img=None,    # ee.Image with inputs as bands (takes precedence o
             NDVI = img['NDVI']
             Tr = img['radiometric_temperature'] # in K
             Ta = img['air_temperature']   # in K
+            Td = img['dewpoint_temperature']   # in K
             P = img['surface_pressure']   # in Pa
             U = img['wind_speed']         # in m/s
             Sdn = img['solar_radiation']  # in W/m2
@@ -226,6 +228,7 @@ def tseb_series(img=None,    # ee.Image with inputs as bands (takes precedence o
             NDVI = to_ndarray(NDVI)
             P = to_ndarray(P)
             Ta = to_ndarray(Ta)
+            Td = to_ndarray(Td)
             U = to_ndarray(U)
             Sdn = to_ndarray(Sdn)
             Ldn = to_ndarray(Ldn)
@@ -255,7 +258,7 @@ def tseb_series(img=None,    # ee.Image with inputs as bands (takes precedence o
                                                                       #                                or exp(-kLAI/sqrt(2cos(Zenith))))
                                                                       # here in TSEB we use the second parameterization (use_zenith = True)
     G = compute_g(doy = doy, time=time, Rns = Rns, G_params = G_params, longitude=longitude) # Soil heat flux (Santanello et al., 2003)
-    met_params = compute_met_params(Ta, P) # computes the following meteorological parameters:
+    met_params = compute_met_params(Ta, Td, P) # computes the following meteorological parameters:
                                              # q: specific humidity (dimensionless)
                                              # ea: water vapor pressure, in Pa
                                              # rho: air density, in kg m-3

tests/test_geeet.py

@@ -21,6 +21,7 @@ def setUp(self):
             NDVI = [0.8, 0.8],
             P = [95500, 95500],
             Ta = [293, 293], 
+            Td = [274, 274],
             U = [5,5], 
             Sdn = [800, 400], 
             Ldn = [300, 200]
@@ -47,6 +48,7 @@ def setUp(self):
                 xr.DataArray(self.tseb_series_inputs["NDVI"]).rename("NDVI"),
                 xr.DataArray(self.tseb_series_inputs["Tr"]).rename("radiometric_temperature"),
                 xr.DataArray(self.tseb_series_inputs["Ta"]).rename("air_temperature"),
+                xr.DataArray(self.tseb_series_inputs["Td"]).rename("dewpoint_temperature"),
                 xr.DataArray(self.tseb_series_inputs["P"]).rename("surface_pressure"),
                 xr.DataArray(self.tseb_series_inputs["U"]).rename("wind_speed"),
                 xr.DataArray(self.tseb_series_inputs["Sdn"]).rename("solar_radiation"),
@@ -61,13 +63,13 @@ def tearDown(self):
 
     def test_tseb_series_list_consistency(self):
         """Test tseb_series outputs consistency (list)."""
-        self.assertCountEqual(self.tseb_series_list["LE"] -
+        self.assertListAlmostEqual(self.tseb_series_list["LE"] -
             self.tseb_series_list["LEs"] - self.tseb_series_list["LEc"],
-            [0,0]
+            [0,0], 10
         )
-        self.assertCountEqual(self.tseb_series_list["Rn"]-
+        self.assertListAlmostEqual(self.tseb_series_list["Rn"]-
             self.tseb_series_list["Rns"] - self.tseb_series_list["Rnc"],
-            [0,0]
+            [0,0], 10
         )
 
     def assertListAlmostEqual(self, list1, list2, tol):
@@ -89,13 +91,13 @@ def test_tseb_series_list_energy_balance(self):
 
     def test_tseb_series_np_consistency(self):
         """Test tseb_series outputs consistency (np.array)."""
-        self.assertCountEqual(self.tseb_series_ndarray["LE"] -
+        self.assertListAlmostEqual(self.tseb_series_ndarray["LE"] -
             self.tseb_series_ndarray["LEs"] - self.tseb_series_ndarray["LEc"],
-            [0,0]
+            [0,0], 10
         )
-        self.assertCountEqual(self.tseb_series_ndarray["Rn"] -
+        self.assertListAlmostEqual(self.tseb_series_ndarray["Rn"] -
             self.tseb_series_ndarray["Rns"] - self.tseb_series_ndarray["Rnc"],
-            [0,0]
+            [0,0],10
         )
 
     def test_tseb_series_np_energy_balance(self):