From 3542bb1b4ce097611f1b5a92d71a1eb5baaa5cb4 Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Sat, 12 Dec 2020 09:16:56 -0800 Subject: [PATCH 1/3] Initial code to support ERA5-Land I still need to check the solar radiation unit conversion, but everything else seemed pretty close. --- README.rst | 36 ++++++++++++ openet/refetgee/__init__.py | 2 +- openet/refetgee/daily.py | 82 ++++++++++++++++++++++++++++ openet/refetgee/hourly.py | 79 ++++++++++++++++++++++++++- openet/refetgee/tests/test_daily.py | 55 +++++++++++++++++++ openet/refetgee/tests/test_hourly.py | 41 ++++++++++++++ 6 files changed, 293 insertions(+), 2 deletions(-) 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 e919068..2decd5b 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.4.0' +__version__ = '0.5.0' diff --git a/openet/refetgee/daily.py b/openet/refetgee/daily.py index cf676fd..979ce39 100644 --- a/openet/refetgee/daily.py +++ b/openet/refetgee/daily.py @@ -684,3 +684,85 @@ def rtma(cls, input_coll, rs=None, zw=None, elev=None, lat=None, method=method, 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 RTMA elevation image + (projects/eddi-noaa/era5_land/elevation) will be used if not set. + lat : ee.Image or ee.Number + Latitude image [degrees]. The latitude will be computed + dynamically using ee.Image.pixelLonLat() 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 W 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/earthengine-legacy/assets/' + 'projects/eddi-noaa/era5_land/elevation')\ + .rename(['elevation']) + if lat is None: + lat = ee.Image('projects/earthengine-legacy/assets/' + 'projects/eddi-noaa/era5_land/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()) + + # TODO: Double check that this is correct + ea_img = calcs._sat_vapor_pressure( + input_coll.select(['dewpoint_temperature_2m']).mean().subtract(273.15)), + # ea_img = calcs._actual_vapor_pressure( + # pair=calcs._air_pressure(elev, method), + # q=input_coll.select(['SPFH']).mean()), + + return cls( + tmax=input_coll.select(['temperature_2m']).max().subtract(273.15), + tmin=input_coll.select(['temperature_2m']).min().subtract(273.15), + ea=ea_img, + # TODO: Check aggregation and units for solar + rs=input_coll.select(['surface_solar_radiation_downwards']).mean() + .divide(3600).multiply(0.0864), + # rs=input_coll.select(['surface_solar_radiation_downwards']).sum() + # .divide(3600 * 24).multiply(0.0864), + 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, + ) \ No newline at end of file diff --git a/openet/refetgee/hourly.py b/openet/refetgee/hourly.py index e39b075..7e31dda 100644 --- a/openet/refetgee/hourly.py +++ b/openet/refetgee/hourly.py @@ -396,4 +396,81 @@ def rtma(cls, input_img, rs=None, zw=None, elev=None, lat=None, lon=None, # time=ee.Number(start_date.getRelative('hour', 'day')), time=ee.Number(start_date.get('hour')), method=method, - ) \ No newline at end of file + ) + + @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 ERA5-Land elevation image + (projects/eddi-noaa/era5_land/elevation) will be used if not set. + lat : ee.Image or ee.Number + Latitude image [degrees]. The latitude will be computed + dynamically using ee.Image.pixelLonLat() if not set. + lon : ee.Image or ee.Number + Longitude image [degrees]. The longitude will be computed + dynamically using ee.Image.pixelLonLat() 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 W 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/earthengine-legacy/assets/' + 'projects/eddi-noaa/era5_land/elevation')\ + .rename(['elevation']) + if lat is None: + lat = ee.Image('projects/earthengine-legacy/assets/' + 'projects/eddi-noaa/era5_land/elevation')\ + .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ + .rename(['latitude']) + if lon is None: + lon = ee.Image('projects/earthengine-legacy/assets/' + 'projects/eddi-noaa/era5_land/elevation')\ + .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ + .rename(['longitude']) + + # TODO: Double check that this is correct + ea_img = calcs._sat_vapor_pressure( + input_img.select(['dewpoint_temperature_2m']).subtract(273.15)), + # ea_img = calcs._actual_vapor_pressure( + # q=input_img.select(['SPFH']), + # pair=calcs._air_pressure(elev, method)), + + return cls( + tmean=input_img.select(['temperature_2m']).subtract(273.15), + ea=ea_img, + # TODO: Check aggregation and units for solar + rs=input_img.select(['surface_solar_radiation_downwards'])\ + .divide(3600).multiply(0.0036), + 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 a0c699b..b164ee2 100644 --- a/openet/refetgee/tests/test_daily.py +++ b/openet/refetgee/tests/test_daily.py @@ -66,9 +66,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']) @@ -81,9 +83,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']) @@ -96,9 +100,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']) @@ -111,9 +117,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']) @@ -127,9 +135,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']) @@ -154,10 +164,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']) @@ -178,10 +190,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) @@ -192,13 +206,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']) @@ -210,13 +227,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']) @@ -230,6 +250,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]) \ @@ -244,9 +265,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']) @@ -311,9 +334,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']) @@ -335,6 +360,36 @@ 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', + 'u_component_of_wind_10m', 'v_component_of_wind_10m'] + + wind_u = d_args['uz'] / (2 ** 0.5) + + era5_coll = ee.ImageCollection.fromImages([ + ee.Image.constant([d_args['tmin'] + 273.15, d_args['tdew'] + 273.15, + d_args['rs'], 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'], d_args['q_asce'], d_args['uz']]) \ + .double().rename(band_names) \ + .set({'system:time_start': ee.Date('2015-07-01T12:00:00', 'UTC').millis()}) + ]) + + refet = Daily.era5( + 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..52d78cd 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'], + wind_u, wind_u])\ + .rename(['temperature_2m', 'dewpoint_temperature_2m', + 'surface_solar_radiation_downwards', + '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']) From 0a6b60ee3d8d714007ec49244e9dae8280aeb5d3 Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Sun, 13 Dec 2020 13:30:58 -0800 Subject: [PATCH 2/3] Fixed the Rs units and calcs --- openet/refetgee/daily.py | 18 ++++++++---------- openet/refetgee/hourly.py | 17 ++++++++--------- openet/refetgee/tests/test_daily.py | 10 ++++++---- openet/refetgee/tests/test_hourly.py | 4 ++-- 4 files changed, 24 insertions(+), 25 deletions(-) diff --git a/openet/refetgee/daily.py b/openet/refetgee/daily.py index 979ce39..0a17c0d 100644 --- a/openet/refetgee/daily.py +++ b/openet/refetgee/daily.py @@ -716,7 +716,7 @@ def era5_land(cls, input_coll, zw=None, elev=None, lat=None, Notes ----- Temperatures are converted from K to C. - Solar radiation is converted from J m-2 to W m-2 to MJ m-2 day-1. + Solar radiation is converted from J m-2 to MJ m-2 day-1. Actual vapor pressure is computed from dew point temperature. """ @@ -727,11 +727,11 @@ def era5_land(cls, input_coll, zw=None, elev=None, lat=None, zw = ee.Number(10) if elev is None: elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/eddi-noaa/era5_land/elevation')\ + 'projects/climate-engine/era5-land/elevation')\ .rename(['elevation']) if lat is None: lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/eddi-noaa/era5_land/elevation')\ + 'projects/climate-engine/era5-land/elevation')\ .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ .rename(['latitude']) @@ -744,20 +744,18 @@ def wind_magnitude(input_img): # TODO: Double check that this is correct ea_img = calcs._sat_vapor_pressure( - input_coll.select(['dewpoint_temperature_2m']).mean().subtract(273.15)), + input_coll.select(['dewpoint_temperature_2m']).mean().subtract(273.15)) # ea_img = calcs._actual_vapor_pressure( # pair=calcs._air_pressure(elev, method), - # q=input_coll.select(['SPFH']).mean()), + # q=input_coll.select(['SPFH']).mean()) return cls( tmax=input_coll.select(['temperature_2m']).max().subtract(273.15), tmin=input_coll.select(['temperature_2m']).min().subtract(273.15), ea=ea_img, - # TODO: Check aggregation and units for solar - rs=input_coll.select(['surface_solar_radiation_downwards']).mean() - .divide(3600).multiply(0.0864), - # rs=input_coll.select(['surface_solar_radiation_downwards']).sum() - # .divide(3600 * 24).multiply(0.0864), + # 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, diff --git a/openet/refetgee/hourly.py b/openet/refetgee/hourly.py index 7e31dda..2c3b491 100644 --- a/openet/refetgee/hourly.py +++ b/openet/refetgee/hourly.py @@ -426,7 +426,7 @@ def era5_land(cls, input_img, zw=None, elev=None, lat=None, lon=None, Notes ----- Temperatures are converted from K to C. - Solar radiation is converted from W m-2 to MJ m-2 day-1. + Solar radiation is converted from J m-2 to MJ m-2 day-1. Actual vapor pressure is computed from dew point temperature. """ @@ -436,32 +436,31 @@ def era5_land(cls, input_img, zw=None, elev=None, lat=None, lon=None, zw = ee.Number(10) if elev is None: elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/eddi-noaa/era5_land/elevation')\ + 'projects/climate-engine/era5-land/elevation')\ .rename(['elevation']) if lat is None: lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/eddi-noaa/era5_land/elevation')\ + 'projects/climate-engine/era5-land/elevation')\ .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ .rename(['latitude']) if lon is None: lon = ee.Image('projects/earthengine-legacy/assets/' - 'projects/eddi-noaa/era5_land/elevation')\ + 'projects/climate-engine/era5-land/elevation')\ .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ .rename(['longitude']) # TODO: Double check that this is correct ea_img = calcs._sat_vapor_pressure( - input_img.select(['dewpoint_temperature_2m']).subtract(273.15)), + input_img.select(['dewpoint_temperature_2m']).subtract(273.15)) # ea_img = calcs._actual_vapor_pressure( # q=input_img.select(['SPFH']), - # pair=calcs._air_pressure(elev, method)), + # pair=calcs._air_pressure(elev, method)) return cls( tmean=input_img.select(['temperature_2m']).subtract(273.15), ea=ea_img, - # TODO: Check aggregation and units for solar - rs=input_img.select(['surface_solar_radiation_downwards'])\ - .divide(3600).multiply(0.0036), + 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']), diff --git a/openet/refetgee/tests/test_daily.py b/openet/refetgee/tests/test_daily.py index b164ee2..4a19060 100644 --- a/openet/refetgee/tests/test_daily.py +++ b/openet/refetgee/tests/test_daily.py @@ -371,22 +371,24 @@ def test_refet_daily_rtma_etr(): 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', + '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, - d_args['rs'], wind_u, wind_u]) \ + 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'], d_args['q_asce'], d_args['uz']]) \ + 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( + 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') diff --git a/openet/refetgee/tests/test_hourly.py b/openet/refetgee/tests/test_hourly.py index 52d78cd..32906bc 100644 --- a/openet/refetgee/tests/test_hourly.py +++ b/openet/refetgee/tests/test_hourly.py @@ -231,10 +231,10 @@ def test_refet_hourly_era5_land_etr(): era5_img = ee.Image.constant([ h_args['tmean'] + 273.15, h_args['tdew'] + 273.15, - h_args['rs'], + h_args['rs'] * 1000000, wind_u, wind_u])\ .rename(['temperature_2m', 'dewpoint_temperature_2m', - 'surface_solar_radiation_downwards', + 'surface_solar_radiation_downwards_hourly', 'u_component_of_wind_10m', 'v_component_of_wind_10m'])\ .set('system:time_start', era5_time) From 30e3a012f2b057c3fb1568d1178b2cac40d63357 Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Mon, 14 Nov 2022 20:02:56 -0800 Subject: [PATCH 3/3] Switching era5land ancillary assets over to openet versions --- .gitignore | 1 + openet/refetgee/daily.py | 35 ++++++++++++--------------- openet/refetgee/hourly.py | 51 ++++++++++++++++++--------------------- 3 files changed, 41 insertions(+), 46 deletions(-) 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/openet/refetgee/daily.py b/openet/refetgee/daily.py index 3a448fe..ef076d5 100644 --- a/openet/refetgee/daily.py +++ b/openet/refetgee/daily.py @@ -783,11 +783,13 @@ def era5_land(cls, input_coll, zw=None, elev=None, lat=None, zw : ee.Number, optional Wind speed height [m] (the default is 10). elev : ee.Image or ee.Number, optional - Elevation image [m]. The RTMA elevation image - (projects/eddi-noaa/era5_land/elevation) will be used if not set. + 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 latitude will be computed - dynamically using ee.Image.pixelLonLat() if not set. + 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. @@ -811,35 +813,30 @@ def era5_land(cls, input_coll, zw=None, elev=None, lat=None, if zw is None: zw = ee.Number(10) if elev is None: - elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/era5-land/elevation')\ + elev = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ .rename(['elevation']) if lat is None: - lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/era5-land/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - .rename(['latitude']) + 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()) - # TODO: Double check that this is correct - ea_img = calcs._sat_vapor_pressure( - input_coll.select(['dewpoint_temperature_2m']).mean().subtract(273.15)) - # ea_img = calcs._actual_vapor_pressure( - # pair=calcs._air_pressure(elev, method), - # q=input_coll.select(['SPFH']).mean()) + 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=ea_img, + 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'])\ + rs=input_coll.select(['surface_solar_radiation_downwards_hourly']) .sum().divide(1000000), uz=wind_img, zw=zw, diff --git a/openet/refetgee/hourly.py b/openet/refetgee/hourly.py index 9cbc865..cad42aa 100644 --- a/openet/refetgee/hourly.py +++ b/openet/refetgee/hourly.py @@ -413,14 +413,17 @@ def era5_land(cls, input_img, zw=None, elev=None, lat=None, lon=None, zw : ee.Number, optional Wind speed height [m] (the default is 10). elev : ee.Image or ee.Number, optional - Elevation image [m]. The ERA5-Land elevation image - (projects/eddi-noaa/era5_land/elevation) will be used if not set. + 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 latitude will be computed - dynamically using ee.Image.pixelLonLat() if not set. + 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 longitude will be computed - dynamically using ee.Image.pixelLonLat() if not set. + 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. @@ -438,34 +441,28 @@ def era5_land(cls, input_img, zw=None, elev=None, lat=None, lon=None, if zw is None: zw = ee.Number(10) if elev is None: - elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/era5-land/elevation')\ + elev = ee.Image('projects/openet/assets/meteorology/era5land/elevation')\ .rename(['elevation']) if lat is None: - lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/era5-land/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - .rename(['latitude']) + 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/earthengine-legacy/assets/' - 'projects/climate-engine/era5-land/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ - .rename(['longitude']) - - # TODO: Double check that this is correct - ea_img = calcs._sat_vapor_pressure( - input_img.select(['dewpoint_temperature_2m']).subtract(273.15)) - # ea_img = calcs._actual_vapor_pressure( - # q=input_img.select(['SPFH']), - # pair=calcs._air_pressure(elev, method)) + 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=ea_img, - rs=input_img.select(['surface_solar_radiation_downwards_hourly'])\ + 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))\ + 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,