diff --git a/.gitignore b/.gitignore index f27e4f7..bffa93e 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ __pycache__/ *$py.class # Distribution / packaging +build/ dist/ .eggs/ *.egg-info/ diff --git a/README.rst b/README.rst index 720446e..a1eca2f 100644 --- a/README.rst +++ b/README.rst @@ -177,6 +177,42 @@ For the daily function, the RTMA collection must be filtered to a single 24 hour print('ETr: {:.2f} mm'.format(float(etr['etr']))) +ERA5-Land +--------- + +Helper functions for computing daily/hourly ETo/ETr for `ERA5-Land `__ images are available. + +For the daily function, the ERA5-Land collection must be filtered to a single 24 hour period. + +.. code-block:: console + + import ee + import openet.refetgee + + source_coll = ee.ImageCollection('ECMWF/ERA5_LAND/HOURLY')\ + .filterDate('2015-07-01', '2015-07-02') + etr = openet.refetgee.Daily.era5_land(source_coll).etr\ + .reduceRegion(reducer=ee.Reducer.first(), + geometry=ee.Geometry.Point(-118.77388, 39.4575), + scale=1000)\ + .getInfo() + + print('ETr: {:.2f} mm'.format(float(etr['etr']))) + +.. code-block:: console + + import ee + import openet.refetgee + + source_img = ee.Image('ECMWF/ERA5_LAND/HOURLY/20150701T20') + etr = openet.refetgee.Hourly.era5_land(source_img).etr\ + .reduceRegion(reducer=ee.Reducer.first(), + geometry=ee.Geometry.Point(-118.77388, 39.4575), + scale=1000)\ + .getInfo() + + print('ETr: {:.2f} mm'.format(float(etr['etr']))) + Input Parameters ================ diff --git a/openet/refetgee/__init__.py b/openet/refetgee/__init__.py index e15a4c8..1b5f283 100644 --- a/openet/refetgee/__init__.py +++ b/openet/refetgee/__init__.py @@ -1,4 +1,4 @@ from .daily import Daily from .hourly import Hourly -__version__ = '0.5.2' +__version__ = '0.6.0' diff --git a/openet/refetgee/daily.py b/openet/refetgee/daily.py index ee6fe37..ef076d5 100644 --- a/openet/refetgee/daily.py +++ b/openet/refetgee/daily.py @@ -770,6 +770,83 @@ def rtma(cls, input_coll, rs=None, zw=None, elev=None, lat=None, rso_type=rso_type, ) + @classmethod + def era5_land(cls, input_coll, zw=None, elev=None, lat=None, + method='asce', rso_type=None): + """Initialize daily RefET from a hourly ERA5-Land image collection + + Parameters + ---------- + input_coll : ee.ImageCollection + Collection of ERA5-Land hourly images for a single day from the + collection ECMWF/ERA5_LAND/HOURLY. + zw : ee.Number, optional + Wind speed height [m] (the default is 10). + elev : ee.Image or ee.Number, optional + Elevation image [m]. The OpenET ERA5-Land elevation image + (projects/openet/assets/meteorology/era5land/elevation) + will be used if not set. + lat : ee.Image or ee.Number + Latitude image [degrees]. The OpenET ERA5-Land latitude image + (projects/openet/assets/meteorology/era5land/latitude) + will be used if not set. + method : {'asce' (default), 'refet'}, optional + Specifies which calculation method to use. + * 'asce' -- Calculations will follow ASCE-EWRI 2005. + * 'refet' -- Calculations will follow RefET software. + rso_type : {None (default), 'full' , 'simple'}, optional + Specifies which clear sky solar radiation (Rso) model to use. + * None -- Rso type will be determined from "method" parameter + * 'full' -- Full clear sky solar formulation + * 'simple' -- Simplified clear sky solar formulation + + Notes + ----- + Temperatures are converted from K to C. + Solar radiation is converted from J m-2 to MJ m-2 day-1. + Actual vapor pressure is computed from dew point temperature. + + """ + input_coll = ee.ImageCollection(input_coll) + start_date = ee.Date(ee.Image(input_coll.first()).get('system:time_start')) + + if zw is None: + zw = ee.Number(10) + if elev is None: + elev = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ + .rename(['elevation']) + if lat is None: + lat = ee.Image('projects/openet/assets/meteorology/era5land/latitude')\ + # lat = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ + # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ + # .rename(['latitude']) + + def wind_magnitude(input_img): + """Compute hourly wind magnitude from vectors""" + return ee.Image(input_img.select(['u_component_of_wind_10m'])).pow(2)\ + .add(ee.Image(input_img.select(['v_component_of_wind_10m'])).pow(2))\ + .sqrt().rename(['wind_10m']) + + wind_img = ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()) + + return cls( + tmax=input_coll.select(['temperature_2m']).max().subtract(273.15), + tmin=input_coll.select(['temperature_2m']).min().subtract(273.15), + ea=calcs._sat_vapor_pressure( + input_coll.select(['dewpoint_temperature_2m']).mean().subtract(273.15) + ), + # TODO: Check that solar does not need additional conversion + rs=input_coll.select(['surface_solar_radiation_downwards_hourly']) + .sum().divide(1000000), + uz=wind_img, + zw=zw, + elev=elev, + lat=lat, + doy=ee.Number(start_date.getRelative('day', 'year')).add(1).double(), + method=method, + rso_type=rso_type, + ) + # @classmethod # def prism(cls, input_img, lat=None): # """Initialize daily RefET from a PRISM image diff --git a/openet/refetgee/hourly.py b/openet/refetgee/hourly.py index 0322b9a..cad42aa 100644 --- a/openet/refetgee/hourly.py +++ b/openet/refetgee/hourly.py @@ -400,3 +400,76 @@ def rtma(cls, input_img, rs=None, zw=None, elev=None, lat=None, lon=None, time=ee.Number(start_date.get('hour')), method=method, ) + + @classmethod + def era5_land(cls, input_img, zw=None, elev=None, lat=None, lon=None, + method='asce'): + """Initialize hourly RefET from an ERA5-Land image + + Parameters + ---------- + input_img : ee.Image + ERA5-Land hourly image from the collection ECMWF/ERA5_LAND/HOURLY. + zw : ee.Number, optional + Wind speed height [m] (the default is 10). + elev : ee.Image or ee.Number, optional + Elevation image [m]. The OpenET ERA5-Land elevation image + (projects/openet/assets/meteorology/era5land/elevation) + will be used if not set. + lat : ee.Image or ee.Number + Latitude image [degrees]. The OpenET ERA5-Land latitude image + (projects/openet/assets/meteorology/era5land/latitude) + will be used if not set. + lon : ee.Image or ee.Number + Longitude image [degrees]. The OpenET ERA5-Land longitude image + (projects/openet/assets/meteorology/era5land/longitude) + will be used if not set. + method : {'asce' (default), 'refet'}, optional + Specifies which calculation method to use. + * 'asce' -- Calculations will follow ASCE-EWRI 2005. + * 'refet' -- Calculations will follow RefET software. + + Notes + ----- + Temperatures are converted from K to C. + Solar radiation is converted from J m-2 to MJ m-2 day-1. + Actual vapor pressure is computed from dew point temperature. + + """ + start_date = ee.Date(input_img.get('system:time_start')) + + if zw is None: + zw = ee.Number(10) + if elev is None: + elev = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ + .rename(['elevation']) + if lat is None: + lat = ee.Image('projects/openet/assets/meteorology/era5land/latitude') + # lat = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ + # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ + # .rename(['latitude']) + if lon is None: + lon = ee.Image('projects/openet/assets/meteorology/era5land/longitude') + # lon = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ + # .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ + # .rename(['longitude']) + + return cls( + tmean=input_img.select(['temperature_2m']).subtract(273.15), + ea=calcs._sat_vapor_pressure( + input_img.select(['dewpoint_temperature_2m']).subtract(273.15) + ), + rs=input_img.select(['surface_solar_radiation_downwards_hourly']) + .divide(1000000), + uz=input_img.select(['u_component_of_wind_10m']).pow(2) + .add(input_img.select(['v_component_of_wind_10m']).pow(2)) + .sqrt().rename(['wind_10m']), + zw=zw, + elev=elev, + lat=lat, + lon=lon, + doy=ee.Number(start_date.getRelative('day', 'year')).add(1).double(), + # time=ee.Number(start_date.getRelative('hour', 'day')), + time=ee.Number(start_date.get('hour')), + method=method, + ) diff --git a/openet/refetgee/tests/test_daily.py b/openet/refetgee/tests/test_daily.py index 2255353..a1cfeff 100644 --- a/openet/refetgee/tests/test_daily.py +++ b/openet/refetgee/tests/test_daily.py @@ -69,9 +69,11 @@ def test_refet_daily_etr(): uz=ee.Image.constant(d_args['uz']), zw=ee.Number(s_args['zw']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='refet') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_refet']) @@ -84,9 +86,11 @@ def test_refet_daily_eto(): uz=ee.Image.constant(d_args['uz']), zw=ee.Number(s_args['zw']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='refet') + output = refet.eto\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['eto']) == pytest.approx(d_args['eto_refet']) @@ -99,9 +103,11 @@ def test_refet_daily_etw(): uz=ee.Image.constant(d_args['uz']), zw=ee.Number(s_args['zw']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='refet') + output = refet.etw\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etw']) == pytest.approx(d_args['etw_refet']) @@ -159,9 +165,11 @@ def test_refet_daily_rso_type_simple(): uz=ee.Image.constant(d_args['uz']), zw=ee.Number(s_args['zw']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='refet', rso_type='simple') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_rso_simple']) @@ -175,9 +183,11 @@ def test_refet_daily_rso_type_array(): elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='refet', rso_type='array', rso=ee.Number(d_args['rso'])) + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_refet']) @@ -202,10 +212,12 @@ def test_refet_daily_etr_asce(): uz=ee.Image.constant(d_args['uz']), zw=ee.Number(s_args['zw']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='asce') + output = refet.etr.reduceRegion( reducer=ee.Reducer.first(), geometry=ee.Geometry.Rectangle([0, 0, 10, 10], 'EPSG:32613', False), scale=1).getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_asce']) @@ -226,10 +238,12 @@ def test_refet_daily_etsz(surface, expected): uz=ee.Image.constant(d_args['uz']), zw=ee.Number(s_args['zw']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), doy=ee.Number(d_args['doy']), method='refet') + output = refet.etsz(surface).rename(['etsz']).reduceRegion( reducer=ee.Reducer.first(), geometry=ee.Geometry.Rectangle([0, 0, 10, 10], 'EPSG:32613', False), scale=1).getInfo() + assert float(output['etsz']) == pytest.approx(expected) @@ -240,13 +254,16 @@ def test_refet_daily_gridmet_etr(): d_args['q_asce'], d_args['rs'] / 0.0864, d_args['uz']])\ .rename(['tmmx', 'tmmn', 'sph', 'srad', 'vs'])\ .set('system:time_start', ee.Date('2015-07-01').millis()) + refet = Daily.gridmet( ee.Image(gridmet_img), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), zw=ee.Number(s_args['zw']), method='asce') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_asce']) @@ -258,13 +275,16 @@ def test_refet_daily_maca_etr(): d_args['uz'] / (2 ** 0.5), d_args['uz'] / (2 ** 0.5)])\ .rename(['tasmax', 'tasmin', 'huss', 'rsds', 'uas', 'vas'])\ .set('system:time_start', ee.Date('2015-07-01').millis()) + refet = Daily.maca( ee.Image(maca_img), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), zw=ee.Number(s_args['zw']), method='asce') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_asce']) @@ -278,6 +298,7 @@ def test_refet_daily_nldas_etr(): 'shortwave_radiation', 'wind_u', 'wind_v'] wind_u = d_args['uz'] / (2 ** 0.5) + nldas_coll = ee.ImageCollection.fromImages([ ee.Image.constant([d_args['tmin'], d_args['q_asce'], 0.0, wind_u, wind_u]) \ @@ -292,9 +313,11 @@ def test_refet_daily_nldas_etr(): refet = Daily.nldas( nldas_coll, elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), zw=ee.Number(s_args['zw']), method='asce') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_asce']) @@ -359,9 +382,11 @@ def test_refet_daily_cfsv2_etr(): refet = Daily.cfsv2( cfsv2_coll, elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), zw=ee.Number(s_args['zw']), method='asce') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(d_args['etr_asce']) @@ -383,6 +408,38 @@ def test_refet_daily_rtma_etr(): rs=ee.Image.constant(d_args['rs']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), zw=ee.Number(s_args['zw']), method='asce') + + output = refet.etr\ + .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ + .getInfo() + + assert float(output['etr']) == pytest.approx(d_args['etr_asce']) + + +def test_refet_daily_era5_land_etr(): + """Generate a fake ERA5-Land image from the test values""" + band_names = ['temperature_2m', 'dewpoint_temperature_2m', + 'surface_solar_radiation_downwards_hourly', + 'u_component_of_wind_10m', 'v_component_of_wind_10m'] + + wind_u = d_args['uz'] / (2 ** 0.5) + + # Each image needs half of the total Rs since hourly images are summed + era5_coll = ee.ImageCollection.fromImages([ + ee.Image.constant([d_args['tmin'] + 273.15, d_args['tdew'] + 273.15, + 0, wind_u, wind_u]) \ + .double().rename(band_names) \ + .set({'system:time_start': ee.Date('2015-07-01T00:00:00', 'UTC').millis()}), + ee.Image.constant([d_args['tmax'] + 273.15, d_args['tdew'] + 273.15, + d_args['rs'] * 1000000, wind_u, wind_u]) \ + .double().rename(band_names) \ + .set({'system:time_start': ee.Date('2015-07-01T12:00:00', 'UTC').millis()}) + ]) + + refet = Daily.era5_land( + era5_coll, + elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), + zw=ee.Number(s_args['zw']), method='asce') output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() diff --git a/openet/refetgee/tests/test_hourly.py b/openet/refetgee/tests/test_hourly.py index 9fd979d..32906bc 100644 --- a/openet/refetgee/tests/test_hourly.py +++ b/openet/refetgee/tests/test_hourly.py @@ -162,10 +162,12 @@ def test_refet_daily_etsz(surface, expected): lat=ee.Number(s_args['lat']), lon=ee.Number(s_args['lon']), doy=ee.Number(h_args['doy']), time=ee.Number(h_args['time']), method='refet') + output = refet.etsz(surface).rename(['etsz']).reduceRegion( reducer=ee.Reducer.first(), geometry=ee.Geometry.Rectangle([0, 0, 10, 10], 'EPSG:32613', False), scale=1).getInfo() + assert float(output['etsz']) == pytest.approx(expected) @@ -173,20 +175,25 @@ def test_refet_hourly_nldas_etr(): """Generate a fake NLDAS image from the test values""" nldas_time = ee.Date( '2015-07-01T{}:00:00'.format(int(h_args['time'])), 'UTC').millis() + wind_u = h_args['uz'] / (2 ** 0.5) + nldas_img = ee.Image.constant([ h_args['tmean'], h_args['q_asce'], h_args['rs'] / 0.0036, wind_u, wind_u])\ .rename(['temperature', 'specific_humidity', 'shortwave_radiation', 'wind_u', 'wind_v'])\ .set('system:time_start', nldas_time) + refet = Hourly.nldas( ee.Image(nldas_img), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), lon=ee.Number(s_args['lon']), zw=ee.Number(s_args['zw']), method='asce') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(h_args['etr_asce']) @@ -194,19 +201,53 @@ def test_refet_hourly_rtma_etr(): """Generate a fake RTMA image from the test values""" rtma_time = ee.Date( '2015-07-01T{}:00:00'.format(int(h_args['time'])), 'UTC').millis() + rtma_img = ee.Image.constant([ h_args['tmean'], h_args['q_asce'], h_args['uz']])\ .rename(['TMP', 'SPFH', 'WIND'])\ .set('system:time_start', rtma_time) + refet = Hourly.rtma( ee.Image(rtma_img), rs=ee.Image.constant(h_args['rs']), elev=ee.Number(s_args['elev']), lat=ee.Number(s_args['lat']), lon=ee.Number(s_args['lon']), zw=ee.Number(s_args['zw']), method='asce') + + output = refet.etr\ + .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ + .getInfo() + + assert float(output['etr']) == pytest.approx(h_args['etr_asce']) + + +def test_refet_hourly_era5_land_etr(): + """Generate a fake ERA5-Land image from the test values""" + era5_time = ee.Date( + '2015-07-01T{}:00:00'.format(int(h_args['time'])), 'UTC').millis() + + wind_u = h_args['uz'] / (2 ** 0.5) + + era5_img = ee.Image.constant([ + h_args['tmean'] + 273.15, + h_args['tdew'] + 273.15, + h_args['rs'] * 1000000, + wind_u, wind_u])\ + .rename(['temperature_2m', 'dewpoint_temperature_2m', + 'surface_solar_radiation_downwards_hourly', + 'u_component_of_wind_10m', 'v_component_of_wind_10m'])\ + .set('system:time_start', era5_time) + + refet = Hourly.era5_land( + ee.Image(era5_img), + elev=ee.Number(s_args['elev']), + lat=ee.Number(s_args['lat']), lon=ee.Number(s_args['lon']), + zw=ee.Number(s_args['zw']), method='asce') + output = refet.etr\ .reduceRegion(ee.Reducer.first(), geometry=constant_geom, scale=1)\ .getInfo() + assert float(output['etr']) == pytest.approx(h_args['etr_asce'])