diff --git a/openet/refetgee/calcs.py b/openet/refetgee/calcs.py index 36fca84..a1a568a 100644 --- a/openet/refetgee/calcs.py +++ b/openet/refetgee/calcs.py @@ -33,8 +33,10 @@ def _air_pressure(elev, method='asce'): if method == 'asce': return elev.multiply(-0.0065).add(293).divide(293).pow(5.26).multiply(101.3) elif method == 'refet': - return elev.multiply(-0.0065).add(293).divide(293)\ + return ( + elev.multiply(-0.0065).add(293).divide(293) .pow(9.8 / (0.0065 * 286.9)).multiply(101.3) + ) def _sat_vapor_pressure(temperature): @@ -55,10 +57,10 @@ def _sat_vapor_pressure(temperature): es = 0.6108 * exp(17.27 * temperature / (temperature + 237.3)) """ - return temperature.add(237.3).pow(-1).multiply(temperature)\ + return ( + temperature.add(237.3).pow(-1).multiply(temperature) .multiply(17.27).exp().multiply(0.6108) - # return temperature.multiply(17.27).divide(temperature.add(237.3)).exp()\ - # .multiply(0.6108) + ) def _es_slope(tmean, method='asce'): @@ -83,11 +85,12 @@ def _es_slope(tmean, method='asce'): """ if method == 'refet': - return _sat_vapor_pressure(tmean)\ - .multiply(4098.0).divide(tmean.add(237.3).pow(2)) + return _sat_vapor_pressure(tmean).multiply(4098.0).divide(tmean.add(237.3).pow(2)) elif method == 'asce': - return tmean.add(237.3).pow(-1).multiply(tmean).multiply(17.27).exp()\ + return ( + tmean.add(237.3).pow(-1).multiply(tmean).multiply(17.27).exp() .multiply(2503.0).divide(tmean.add(237.3).pow(2)) + ) def _actual_vapor_pressure(q, pair): @@ -173,7 +176,6 @@ def _precipitable_water(ea, pair): ee.Image or ee.Number Precipitable water [mm]. - Notes ----- w = pair * 0.14 * ea + 2.1 @@ -213,7 +215,7 @@ def _delta(doy, method='asce'): Returns ------- - ee.Image or ee.Number + delta : ee.Image or ee.Number Earth declination [radians]. Notes @@ -225,10 +227,14 @@ def _delta(doy, method='asce'): """ if method == 'asce': - return _doy_fraction(doy).subtract(1.39).sin().multiply(0.409) + delta = _doy_fraction(doy).subtract(1.39).sin().multiply(0.409) else: - return doy.add(284).multiply(2 * math.pi / 365).sin()\ + delta = ( + doy.add(284).multiply(2 * math.pi / 365).sin() .multiply(23.45 * (math.pi / 180)) + ) + + return delta def _dr(doy): @@ -272,8 +278,10 @@ def _seasonal_correction(doy): """ b = doy.subtract(81).divide(364.0).multiply(2 * math.pi) - return b.multiply(2).sin().multiply(0.1645)\ + return ( + b.multiply(2).sin().multiply(0.1645) .subtract(b.cos().multiply(0.1255)).subtract(b.sin().multiply(0.0250)) + ) def _solar_time_rad(lon, time_mid, sc): @@ -323,8 +331,7 @@ def _omega(solar_time): # Need to adjust omega so that the values go from -pi to pi # Values outside this range are wrapped (i.e. -3*pi/2 -> pi/2) - omega = _wrap(omega, -math.pi, math.pi) - return omega + return _wrap(omega, -math.pi, math.pi) def _wrap(x, x_min, x_max): @@ -353,6 +360,7 @@ def _wrap(x, x_min, x_max): """ x_range = x_max - x_min + return x.subtract(x_min).mod(x_range).add(x_range).mod(x_range).add(x_min) # return np.mod((x - x_min), (x_max - x_min)) + x_min @@ -407,8 +415,10 @@ def _ra_daily(lat, doy, method='asce'): """ delta = _delta(doy, method) omegas = _omega_sunset(lat, delta) - theta = omegas.multiply(lat.sin()).multiply(delta.sin())\ + theta = ( + omegas.multiply(lat.sin()).multiply(delta.sin()) .add(lat.cos().multiply(delta.cos()).multiply(omegas.sin())) + ) if method == 'asce': # (24. / math.pi) * 4.92 * _dr(doy) * theta @@ -416,6 +426,7 @@ def _ra_daily(lat, doy, method='asce'): else: # ra = (24. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta ra = theta.multiply(_dr(doy)).multiply((24. / math.pi) * (1367 * 0.0036)) + return ra @@ -467,6 +478,7 @@ def _ra_hourly(lat, lon, doy, time_mid, method='asce'): else: # ra = (12. / math.pi) * (1367 * 0.0036) * _dr(doy) * theta ra = theta.multiply(_dr(doy)).multiply((12. / math.pi) * (1367 * 0.0036)) + return ra @@ -516,8 +528,7 @@ def _rso_daily(ea, pair, ra, doy, lat): # Transmissivity index for diffuse radiation (Eq. D.4) kd = kb.multiply(-0.36).add(0.35).min(kb.multiply(0.82).add(0.18)) - rso = kb.add(kd).multiply(ra) - return rso + return kb.add(kd).multiply(ra) def _rso_hourly(ea, ra, pair, doy, time_mid, lat, lon, method='asce'): @@ -557,8 +568,10 @@ def _rso_hourly(ea, ra, pair, doy, time_mid, lat, lon, method='asce'): # sin of the angle of the sun above the horizon (D.6 and Eq. 62) delta = _delta(doy, method) - sin_beta = lat.sin().multiply(delta.sin())\ + sin_beta = ( + lat.sin().multiply(delta.sin()) .add(lat.cos().multiply(delta.cos()).multiply(omega.cos())) + ) # Precipitable water w = _precipitable_water(ea, pair) @@ -572,15 +585,14 @@ def _rso_hourly(ea, ra, pair, doy, time_mid, lat, lon, method='asce'): .exp().multiply(0.98) ) # kb = ea.expression( - # '0.98 * exp((-0.00146 * pair) / (kt * sin_beta) - ' - # ' 0.075 * (w / sin_beta) ** 0.4))', - # {'pair': pair, 'kt': kt, 'sin_beta': sin_beta.max(0.01), 'w': w}) + # '0.98 * exp((-0.00146 * pair) / (kt * sin_beta) - 0.075 * (w / sin_beta) ** 0.4))', + # {'pair': pair, 'kt': kt, 'sin_beta': sin_beta.max(0.01), 'w': w} + # ) # Transmissivity index for diffuse radiation (Eq. D.4) kd = kb.multiply(-0.36).add(0.35).min(kb.multiply(0.82).add(0.18)) - rso = kb.add(kd).multiply(ra) - return rso + return kb.add(kd).multiply(ra) def _rso_simple(ra, elev): @@ -663,8 +675,10 @@ def _fcd_hourly(rs, rso, doy, time_mid, lat, lon, method='asce'): sc = _seasonal_correction(doy) delta = _delta(doy, method) omega = _omega(_solar_time_rad(lon, time_mid, sc)) - beta = lat.sin().multiply(delta.sin())\ + beta = ( + lat.sin().multiply(delta.sin()) .add(lat.cos().multiply(delta.cos()).multiply(omega.cos())).asin() + ) fcd = rs.divide(rso).max(0.3).min(1).multiply(1.35).subtract(0.35) diff --git a/openet/refetgee/daily.py b/openet/refetgee/daily.py index 3b33e5c..f3722af 100644 --- a/openet/refetgee/daily.py +++ b/openet/refetgee/daily.py @@ -116,9 +116,11 @@ def __init__(self, tmax, tmin, ea, rs, uz, zw, elev, lat, doy, method='asce', rs self.es_slope = calcs._es_slope(self.tmean, method) # Saturated vapor pressure - self.es = calcs._sat_vapor_pressure(self.tmax)\ - .add(calcs._sat_vapor_pressure(self.tmin))\ + self.es = ( + calcs._sat_vapor_pressure(self.tmax) + .add(calcs._sat_vapor_pressure(self.tmin)) .multiply(0.5) + ) # Vapor pressure deficit self.vpd = calcs._vpd(es=self.es, ea=self.ea) @@ -134,13 +136,13 @@ def __init__(self, tmax, tmin, ea, rs, uz, zw, elev, lat, doy, method='asce', rs self.rso = calcs._rso_simple(ra=self.ra, elev=self.elev) elif method.lower() == 'refet': self.rso = calcs._rso_daily( - ea=self.ea, pair=self.pair, ra=self.ra, doy=self.doy, lat=self.lat + ea=self.ea, pair=self.pair, ra=self.ra, doy=self.doy, lat=self.lat, ) elif rso_type.lower() == 'simple': self.rso = calcs._rso_simple(ra=self.ra, elev=self.elev) elif rso_type.lower() == 'full': self.rso = calcs._rso_daily( - ea=self.ea, pair=self.pair, ra=self.ra, doy=self.doy, lat=self.lat + ea=self.ea, pair=self.pair, ra=self.ra, doy=self.doy, lat=self.lat, ) elif rso_type.lower() == 'array': # Use rso array passed to function @@ -211,9 +213,12 @@ def _etsz(self): # return self.tmin.expression( # '(0.408 * es_slope * rn + (psy * cn * u2 * vpd / (tmean + 273))) / ' # '(es_slope + psy * (cd * u2 + 1))', - # {'cd': self.cd, 'cn': self.cn, 'es_slope': self.es_slope, - # 'psy': self.psy, 'rn': self.rn, 'tmean': self.tmean, - # 'u2': self.u2, 'vpd': self.vpd}) + # { + # 'cd': self.cd, 'cn': self.cn, 'es_slope': self.es_slope, + # 'psy': self.psy, 'rn': self.rn, 'tmean': self.tmean, + # 'u2': self.u2, 'vpd': self.vpd, + # } + # ) @lazy_property def etw(self): @@ -233,15 +238,17 @@ def etw(self): self.es_slope .expression( '(alpha * es_slope * rn * 1000 / (2453 * (es_slope + psy)))', - {'es_slope': self.es_slope, 'rn': self.rn, 'psy': self.psy, 'alpha': 1.26}) + {'es_slope': self.es_slope, 'rn': self.rn, 'psy': self.psy, 'alpha': 1.26} + ) .rename(['etw']) .set('system:time_start', self.time_start) ) - # return self.es_slope.multiply(self.rn)\ - # .divide((self.es_slope.add(self.psy)).multiply(2453))\ - # .multiply(1.26).multiply(1000)\ - # .rename(['etw'])\ + # return ( + # self.es_slope.add(self.psy).multiply(2453).pow(-1).multiply(1.26 * 1000) + # .multiply(self.es_slope).multiply(self.rn).multiply(1.26) + # .rename(['etw']) # .set('system:time_start', self.time_start) + # ) @lazy_property def eto_fs1(self): @@ -261,14 +268,17 @@ def eto_fs1(self): self.u2 .expression( '(delta / (delta + psy * (1 + 0.34 * u2))) * (0.408 * rn)', - {'delta': self.es_slope, 'psy': self.psy, 'u2': self.u2, 'rn': self.rn}) + {'delta': self.es_slope, 'psy': self.psy, 'u2': self.u2, 'rn': self.rn} + ) .rename(['eto_fs1']) .set('system:time_start', self.time_start) ) - - # return self.es_slope\ - # .divide(self.es_slope.add(self.psy.multiply(self.u2.multiply(0.34).add(1))))\ - # .multiply(self.rn.multiply(0.408)) + # return ( + # self.u2.multiply(0.34).add(1).multiply(self.psy).add(self.es_slope) + # .pow(-1).multiply(self.es_slope).multiply(self.rn).multiply(0.408) + # .rename(['eto_fs1']) + # .set('system:time_start', self.time_start) + # ) @lazy_property def eto_fs2(self): @@ -286,6 +296,7 @@ def eto_fs2(self): """ # Temperature Term (Eq. 14) tt = self.u2.expression('(900 / (t + 273)) * u2', {'t': self.tmean, 'u2': self.u2}) + # Psi Term (Eq. 13) pt = self.u2.expression( 'psy / (slope + psy * (1 + 0.34 * u2))', @@ -293,15 +304,10 @@ def eto_fs2(self): ) return ( - self.u2 - .expression('PT * TT * (es-ea)', {'PT': pt, 'TT': tt, 'es': self.es, 'ea': self.ea}) + self.es.subtract(self.ea).multiply(pt).multiply(tt) .rename(['eto_fs2']) .set('system:time_start', self.time_start) ) - # return self.PT.multiply(self.TT)\ - # .multiply(self.es.subtract(self.ea))\ - # .rename(['eto_fs2'])\ - # .set('system:time_start', self.time_start) @lazy_property def pet_hargreaves(self): @@ -320,16 +326,17 @@ def pet_hargreaves(self): self.tmax .expression( '0.0023 * (tmean + 17.8) * ((tmax - tmin) ** 0.5) * 0.408 * ra', - {'tmean': self.tmean, 'tmax': self.tmax, 'tmin': self.tmin, 'ra': self.ra}) + {'tmean': self.tmean, 'tmax': self.tmax, 'tmin': self.tmin, 'ra': self.ra} + ) .rename(['pet_hargreaves']) .set('system:time_start', self.time_start) ) - # return self.ra\ - # .multiply(self.tmean.add(17.8))\ - # .multiply(self.tmax.subtract(self.tmin).pow(0.5))\ - # .multiply(0.0023 * 0.408)\ - # .rename(['pet_hargreaves'])\ + # return ( + # self.tmax.subtract(self.tmin).pow(0.5).multiply().multiply(0.0023 * 0.408) + # .multiply(self.tmean.add(17.8)) + # .rename(['pet_hargreaves']) # .set('system:time_start', self.time_start) + # ) @classmethod def gridmet(cls, input_img, zw=None, elev=None, lat=None, method='asce', rso_type=None): @@ -343,7 +350,7 @@ def gridmet(cls, input_img, zw=None, elev=None, lat=None, method='asce', rso_typ Wind speed height [m] (the default is 10). elev : ee.Image or ee.Number, optional Elevation image [m]. The standard GRIDMET elevation image - (projects/climate-engine/gridmet/elevtion) will be used if not set. + (projects/openet/assets/meteorology/gridmet/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. @@ -367,34 +374,20 @@ def gridmet(cls, input_img, zw=None, elev=None, lat=None, method='asce', rso_typ """ image_date = ee.Date(input_img.get('system:time_start')) - # gridmet_transform = [0.041666666666666664, 0, -124.78749996666667, - # 0, -0.041666666666666664, 49.42083333333334] - if zw is None: zw = ee.Number(10) if elev is None: - elev = ( - ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/gridmet/elevation')\ - .rename(['elevation']) - ) - # elev = ee.Image('CGIAR/SRTM90_V4').reproject('EPSG:4326', gridmet_transform) + elev = ee.Image('projects/openet/assets/meteorology/gridmet/elevation') if lat is None: - lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/gridmet/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - .rename(['latitude']) - # lat = ee.Image.pixelLonLat().select('latitude')\ - # .reproject('EPSG:4326', gridmet_transform) - # lat = input_img.select([0]).multiply(0)\ - # .add(ee.Image.pixelLonLat().select('latitude')) + lat = ee.Image('projects/openet/assets/meteorology/gridmet/latitude') return cls( tmax=input_img.select(['tmmx']).subtract(273.15), tmin=input_img.select(['tmmn']).subtract(273.15), ea=calcs._actual_vapor_pressure( + pair=calcs._air_pressure(elev, method), q=input_img.select(['sph']), - pair=calcs._air_pressure(elev, method)), + ), rs=input_img.select(['srad']).multiply(0.0864), uz=input_img.select(['vs']), zw=zw, @@ -417,7 +410,7 @@ def maca(cls, input_img, zw=None, elev=None, lat=None, method='asce', rso_type=N Wind speed height [m] (the default is 10). elev : ee.Image or ee.Number, optional Elevation image [m]. The standard GRIDMET elevation image - (projects/climate-engine/gridmet/elevtion) will be used if not set. + (projects/openet/assets/meteorology/gridmet/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. @@ -442,23 +435,22 @@ def maca(cls, input_img, zw=None, elev=None, lat=None, method='asce', rso_type=N """ image_date = ee.Date(input_img.get('system:time_start')) - # CHECK MACA transform (Why the difference?) + # TODO: CHECK MACA transform (Why the difference?) # maca_transform = [0.04163593073593073, 0, -124.7722, - # 0, 0.04159470085470086, 25.0631] + # 0, 0.04159470085470086, 25.0631] # gridmet_transform = [0.041666666666666664, 0, -124.78749996666667, # 0, -0.041666666666666664, 49.42083333333334] if zw is None: zw = ee.Number(10) if elev is None: - # should we use gridmet elev? maca grid is shifted - elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/gridmet/elevation') + # TODO: should we use gridmet elev? maca grid is shifted + elev = ee.Image('projects/openet/assets/meteorology/gridmet/elevation') # elev = ee.Image('CGIAR/SRTM90_V4').reproject('EPSG:4326', gridmet_transform) if lat is None: - lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/gridmet/elevation') \ - .multiply(0).add(ee.Image.pixelLonLat().select('latitude')) + elev = ee.Image('projects/openet/assets/meteorology/gridmet/latitude') + # lat = ee.Image('projects/openet/assets/meteorology/gridmet/elevation') \ + # .multiply(0).add(ee.Image.pixelLonLat().select('latitude')) # lat = ee.Image.pixelLonLat().select('latitude') \ # .reproject('EPSG:4326', gridmet_transform) # lat = input_img.select([0]).multiply(0) \ @@ -466,19 +458,21 @@ def maca(cls, input_img, zw=None, elev=None, lat=None, method='asce', rso_type=N def wind_magnitude(input_img): """Compute daily wind magnitude from vectors""" - return ee.Image(input_img.select(['uas'])).pow(2) \ - .add(ee.Image(input_img.select(['vas'])).pow(2)) \ + return ( + ee.Image(input_img.select(['uas'])).pow(2) + .add(ee.Image(input_img.select(['vas'])).pow(2)) .sqrt().rename(['uz']) - wind_img = ee.Image(wind_magnitude(input_img)) + ) return cls( tmax=input_img.select(['tasmax']).subtract(273.15), tmin=input_img.select(['tasmin']).subtract(273.15), ea=calcs._actual_vapor_pressure( pair=calcs._air_pressure(elev, method), - q=input_img.select(['huss'])), + q=input_img.select(['huss']), + ), rs=input_img.select(['rsds']).multiply(0.0864), - uz=wind_img.select(['uz']), + uz=ee.Image(wind_magnitude(input_img)).select(['uz']), zw=zw, elev=elev, lat=lat, @@ -528,37 +522,26 @@ def nldas(cls, input_coll, zw=None, elev=None, lat=None, method='asce', rso_type zw = ee.Number(10) if elev is None: elev = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/elevation') - # elev = ee.Image('CGIAR/SRTM90_V4')\ - # .reproject('EPSG:4326', [0.125, 0, -125, 0, -0.125, 53]) if lat is None: lat = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/latitude') - # lat = ee.Image('projects/earthengine-legacy/assets/' - # 'projects/eddi-noaa/nldas/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - # .rename(['latitude']) - # lat = ee.Image.pixelLonLat().select('latitude')\ - # .reproject('EPSG:4326', [0.125, 0, -125, 0, -0.125, 53]) - # lat = input_coll.first().select([0]).multiply(0)\ - # .add(ee.Image.pixelLonLat().select('latitude')) def wind_magnitude(input_img): """Compute hourly wind magnitude from vectors""" - return ee.Image(input_img.select(['wind_u'])).pow(2)\ - .add(ee.Image(input_img.select(['wind_v'])).pow(2))\ + return ( + ee.Image(input_img.select(['wind_u'])).pow(2) + .add(ee.Image(input_img.select(['wind_v'])).pow(2)) .sqrt().rename(['uz']) - wind_img = ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()) - - ea_img = calcs._actual_vapor_pressure( - pair=calcs._air_pressure(elev, method), - q=input_coll.select(['specific_humidity']).mean() - ) + ) return cls( tmax=input_coll.select(['temperature']).max(), tmin=input_coll.select(['temperature']).min(), - ea=ea_img, + ea=calcs._actual_vapor_pressure( + pair=calcs._air_pressure(elev, method), + q=input_coll.select(['specific_humidity']).mean(), + ), rs=input_coll.select(['shortwave_radiation']).sum().multiply(0.0036), - uz=wind_img, + uz=ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()), zw=zw, elev=elev, lat=lat, @@ -628,25 +611,19 @@ def wind_magnitude(input_img): v_img = ee.Image(input_img).select(['v-component_of_wind_height_above_ground']) return u_img.pow(2).add(v_img.pow(2)).sqrt() - wind_img = ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()) - - ea_img = calcs._actual_vapor_pressure( - pair=calcs._air_pressure(elev, method), - # pair=input_coll.select(['Pressure_surface']).mean() - q=input_coll.select(['Specific_humidity_height_above_ground']).mean() - ) - return cls( tmax=input_coll.select(['Maximum_temperature_height_above_ground_6_Hour_Interval']) .max().subtract(273.15), tmin=input_coll.select(['Minimum_temperature_height_above_ground_6_Hour_Interval']) .min().subtract(273.15), - ea=ea_img, + ea=calcs._actual_vapor_pressure( + pair=calcs._air_pressure(elev, method), + q=input_coll.select(['Specific_humidity_height_above_ground']).mean(), + ), # TODO: Check the conversion on solar - rs=input_coll - .select(['Downward_Short-Wave_Radiation_Flux_surface_6_Hour_Average']) + rs=input_coll.select(['Downward_Short-Wave_Radiation_Flux_surface_6_Hour_Average']) .mean().multiply(0.0864), - uz=wind_img, + uz=ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()), zw=zw, elev=elev, lat=lat, @@ -700,7 +677,7 @@ def rtma(cls, input_coll, rs=None, zw=None, elev=None, lat=None, method='asce', pass elif isinstance(rs, ee.Number) or isinstance(rs, float) or isinstance(rs, int): rs = ee.Image.constant(rs) - elif rs is None or rs.upper() == 'GRIDMET': + elif (rs is None) or (rs.upper() == 'GRIDMET'): rs_coll = ( ee.ImageCollection('IDAHO_EPSCOR/GRIDMET') .filterDate(start_date, start_date.advance(1, 'day')) @@ -726,35 +703,32 @@ def rtma(cls, input_coll, rs=None, zw=None, elev=None, lat=None, method='asce', zw = ee.Number(10) if elev is None: elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/rtma/elevation')\ - .rename(['elevation']) + 'projects/climate-engine/rtma/elevation') if lat is None: - lat = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/rtma/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - .rename(['latitude']) + lat = ( + ee.Image('projects/earthengine-legacy/assets/' + 'projects/climate-engine/rtma/elevation') + .multiply(0).add(ee.Image.pixelLonLat().select('latitude')) + ) - # DEADBEEF - Don't compute wind speed from the components since - # a wind speed band is provided + # # Not computing wind speed from the components since wind speed is available # def wind_magnitude(input_img): # """Compute hourly wind magnitude from vectors""" # return ee.Image(input_img.select(['UGRD'])).pow(2)\ # .add(ee.Image(input_img.select(['VGRD'])).pow(2))\ # .sqrt().rename(['WIND']) - # wind_img = ee.Image( - # ee.ImageCollection(input_coll.map(wind_magnitude)).mean()) - - ea_img = calcs._actual_vapor_pressure( - pair=calcs._air_pressure(elev, method), - q=input_coll.select(['SPFH']).mean() - ) return cls( tmax=input_coll.select(['TMP']).max(), tmin=input_coll.select(['TMP']).min(), - ea=ea_img, + ea=calcs._actual_vapor_pressure( + pair=calcs._air_pressure(elev, method), + q=input_coll.select(['SPFH']).mean(), + ), rs=rs, uz=input_coll.select(['WIND']).mean(), + # Not computing wind speed from the components since wind speed is available + # uz=ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()) zw=zw, elev=elev, lat=lat, @@ -807,10 +781,7 @@ def era5(cls, input_coll, zw=None, elev=None, lat=None, method='asce', rso_type= if elev is None: elev = ee.Image('projects/openet/assets/meteorology/era5/ancillary/elevation') if lat is None: - lat = ee.Image('projects/openet/assets/meteorology/era5/ancillary/latitude')\ - # lat = ee.Image('projects/openet/assets/meteorology/era5/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - # .rename(['latitude']) + lat = ee.Image('projects/openet/assets/meteorology/era5/ancillary/latitude') def wind_magnitude(input_img): """Compute hourly wind magnitude from vectors""" @@ -818,17 +789,14 @@ def wind_magnitude(input_img): .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) ), - rs=input_coll.select(['surface_solar_radiation_downwards']) - .sum().divide(1000000), - uz=wind_img, + rs=input_coll.select(['surface_solar_radiation_downwards']).sum().divide(1000000), + uz=ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()), zw=zw, elev=elev, lat=lat, @@ -881,18 +849,15 @@ def era5_land(cls, input_coll, zw=None, elev=None, lat=None, method='asce', rso_ if elev is None: elev = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/elevation') if lat is None: - lat = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/latitude')\ - # lat = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - # .rename(['latitude']) + lat = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/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))\ + 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), @@ -900,9 +865,8 @@ def wind_magnitude(input_img): ea=calcs._sat_vapor_pressure( input_coll.select(['dewpoint_temperature_2m']).mean().subtract(273.15) ), - rs=input_coll.select(['surface_solar_radiation_downwards_hourly']) - .sum().divide(1000000), - uz=wind_img, + rs=input_coll.select(['surface_solar_radiation_downwards_hourly']).sum().divide(1000000), + uz=ee.Image(ee.ImageCollection(input_coll.map(wind_magnitude)).mean()), zw=zw, elev=elev, lat=lat, diff --git a/openet/refetgee/hourly.py b/openet/refetgee/hourly.py index dbb1dde..faa8f96 100644 --- a/openet/refetgee/hourly.py +++ b/openet/refetgee/hourly.py @@ -112,7 +112,7 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, method='asc # Extraterrestrial radiation time_mid = self.time.add(0.5) self.ra = calcs._ra_hourly( - lat=self.lat, lon=self.lon, doy=self.doy, time_mid=time_mid, method=method + lat=self.lat, lon=self.lon, doy=self.doy, time_mid=time_mid, method=method, ) # Clear sky solar radiation @@ -121,7 +121,7 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, method='asc elif method == 'refet': self.rso = calcs._rso_hourly( ea=self.ea, ra=self.ra, pair=self.pair, doy=self.doy, - time_mid=time_mid, lat=self.lat, lon=self.lon, method=method + time_mid=time_mid, lat=self.lat, lon=self.lon, method=method, ) # Cloudiness fraction @@ -131,7 +131,7 @@ def __init__(self, tmean, ea, rs, uz, zw, elev, lat, lon, doy, time, method='asc # Beta (not SinBeta) is used for clamping fcd. self.fcd = calcs._fcd_hourly( rs=self.rs, rso=self.rso, doy=self.doy, time_mid=self.time, - lat=self.lat, lon=self.lon, method=method + lat=self.lat, lon=self.lon, method=method, ) # Net long-wave radiation @@ -214,11 +214,6 @@ def _etsz(self): Standardized reference ET [mm]. """ - # return self.u2.multiply(self.vpd).multiply(self.cn).multiply(self.psy)\ - # .divide(self.tmean.add(273))\ - # .add(self.es_slope.multiply(self.rn.subtract(self.g)).multiply(0.408))\ - # .divide(self.u2.multiply(self.cd).add(1)\ - # .multiply(self.psy).add(self.es_slope)) return self.tmean.expression( '(0.408 * es_slope * (rn - g) + (psy * cn * u2 * vpd / (tmean + 273))) / ' @@ -226,9 +221,15 @@ def _etsz(self): { 'cd': self.cd, 'cn': self.cn, 'es_slope': self.es_slope, 'g': self.g, 'psy': self.psy, 'rn': self.rn, 'tmean': self.tmean, - 'u2': self.u2, 'vpd': self.vpd + 'u2': self.u2, 'vpd': self.vpd, } ) + # return ( + # self.u2.multiply(self.vpd).multiply(self.cn).multiply(self.psy) + # .divide(self.tmean.add(273)) + # .add(self.es_slope.multiply(self.rn.subtract(self.g)).multiply(0.408)) + # .divide(self.u2.multiply(self.cd).add(1).multiply(self.psy).add(self.es_slope)) + # ) @classmethod def nldas(cls, input_img, zw=None, elev=None, lat=None, lon=None, method='asce'): @@ -270,32 +271,17 @@ def nldas(cls, input_img, zw=None, elev=None, lat=None, lon=None, method='asce') zw = ee.Number(10) if elev is None: elev = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/elevation') - # elev = ee.Image('CGIAR/SRTM90_V4')\ - # .reproject('EPSG:4326', [0.125, 0, -125, 0, -0.125, 53]) if lat is None: lat = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/latitude') - # lat = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - # .rename(['latitude']) - # lat = ee.Image.pixelLonLat().select('latitude')\ - # .reproject('EPSG:4326', [0.125, 0, -125, 0, -0.125, 53]) - # lat = nldas_img.select([0]).multiply(0)\ - # .add(ee.Image.pixelLonLat().select('latitude')) if lon is None: lon = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/longitude') - # lon = ee.Image('projects/openet/assets/meteorology/nldas/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ - # .rename(['longitude']) - # lon = ee.Image.pixelLonLat().select('longitude')\ - # .reproject('EPSG:4326', [0.125, 0, -125, 0, -0.125, 53]) - # lon = nldas_img.select([0]).multiply(0)\ - # .add(ee.Image.pixelLonLat().select('longitude')) return cls( tmean=input_img.select(['temperature']), ea=calcs._actual_vapor_pressure( + pair=calcs._air_pressure(elev, method), q=input_img.select(['specific_humidity']), - pair=calcs._air_pressure(elev, method)), + ), rs=input_img.select(['shortwave_radiation']).multiply(0.0036), uz=input_img.select(['wind_u']).pow(2).add(input_img.select(['wind_v']).pow(2)) .sqrt().rename(['uz']), @@ -350,10 +336,12 @@ def rtma(cls, input_img, rs=None, zw=None, elev=None, lat=None, lon=None, method pass elif isinstance(rs, ee.Number) or isinstance(rs, float) or isinstance(rs, int): rs = ee.Image.constant(rs) - elif rs is None or rs.upper() == 'NLDAS': - rs = ee.ImageCollection('NASA/NLDAS/FORA0125_H002')\ - .filterDate(start_date, start_date.advance(30, 'minute'))\ + elif (rs is None) or (rs.upper() == 'NLDAS'): + rs = ( + ee.ImageCollection('NASA/NLDAS/FORA0125_H002') + .filterDate(start_date, start_date.advance(30, 'minute')) .select(['shortwave_radiation']) + ) rs = ee.Image(rs.first()).multiply(0.0036) else: raise ValueError('Unsupported Rs input') @@ -362,30 +350,25 @@ def rtma(cls, input_img, rs=None, zw=None, elev=None, lat=None, lon=None, method zw = ee.Number(10) if elev is None: elev = ee.Image('projects/earthengine-legacy/assets/' - 'projects/climate-engine/rtma/elevation')\ - .rename(['elevation']) + 'projects/climate-engine/rtma/elevation') if lat is None: lat = ee.Image('projects/earthengine-legacy/assets/' 'projects/climate-engine/rtma/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - .rename(['latitude']) + .multiply(0).add(ee.Image.pixelLonLat().select('latitude')) if lon is None: lon = ee.Image('projects/earthengine-legacy/assets/' 'projects/climate-engine/rtma/elevation')\ - .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ - .rename(['longitude']) + .multiply(0).add(ee.Image.pixelLonLat().select('longitude')) return cls( tmean=input_img.select(['TMP']), ea=calcs._actual_vapor_pressure( - q=input_img.select(['SPFH']), - pair=calcs._air_pressure(elev, method)), + pair=calcs._air_pressure(elev, method), q=input_img.select(['SPFH']), + ), rs=rs, # Use wind speed band directly instead of computing from components uz=input_img.select(['WIND']), - # uz=input_img.select(['UGRD']).pow(2)\ - # .add(input_img.select(['VGRD']).pow(2))\ - # .sqrt().rename(['WIND']), + # uz=input_img.select(['UGRD']).pow(2).add(input_img.select(['VGRD']).pow(2)).sqrt(), zw=zw, elev=elev, lat=lat, @@ -438,14 +421,8 @@ def era5(cls, input_img, zw=None, elev=None, lat=None, lon=None, method='asce'): elev = ee.Image('projects/openet/assets/meteorology/era5/ancillary/elevation') if lat is None: lat = ee.Image('projects/openet/assets/meteorology/era5/ancillary/latitude') - # lat = ee.Image('projects/openet/assets/meteorology/era5/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - # .rename(['latitude']) if lon is None: lon = ee.Image('projects/openet/assets/meteorology/era5/ancillary/longitude') - # lon = ee.Image('projects/openet/assets/meteorology/era5/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('longitude'))\ - # .rename(['longitude']) return cls( tmean=input_img.select(['temperature_2m']).subtract(273.15), @@ -508,22 +485,15 @@ def era5_land(cls, input_img, zw=None, elev=None, lat=None, lon=None, method='as elev = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/elevation') if lat is None: lat = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/latitude') - # lat = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/elevation')\ - # .multiply(0).add(ee.Image.pixelLonLat().select('latitude'))\ - # .rename(['latitude']) if lon is None: lon = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/longitude') - # lon = ee.Image('projects/openet/assets/meteorology/era5land/ancillary/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), + 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_calcs.py b/openet/refetgee/tests/test_calcs.py index b003010..9db199c 100644 --- a/openet/refetgee/tests/test_calcs.py +++ b/openet/refetgee/tests/test_calcs.py @@ -413,8 +413,8 @@ def test_ra_hourly_refet(lat=s_args['lat'], lon=s_args['lon'], assert float(output) == pytest.approx(expected) def test_ra_hourly_position(lat=s_args['lat'], lon=s_args['lon'], - doy=h_args['doy'], time_mid=h_args['time_mid'], - expected=h_args['ra_asce']): + doy=h_args['doy'], time_mid=h_args['time_mid'], + expected=h_args['ra_asce']): output = utils.get_info(calcs._ra_hourly( ee.Number(lat), ee.Number(lon), ee.Number(doy), ee.Number(time_mid) )) @@ -621,7 +621,7 @@ def test_fcd_hourly_position(rs=h_args['rs'], rso=h_args['rso'], expected=h_args['fcd_asce']): output = utils.get_info(calcs._fcd_hourly( ee.Number(rs), ee.Number(rso), ee.Number(doy), ee.Number(time_mid), - ee.Number(lat), ee.Number(lon) + ee.Number(lat), ee.Number(lon), )) assert float(output) == pytest.approx(expected) @@ -641,7 +641,7 @@ def test_rnl_daily_number(tmin=d_args['tmin'], tmax=d_args['tmax'], ea=d_args['ea'], fcd=d_args['fcd'], expected=d_args['rnl']): output = utils.get_info(calcs._rnl_daily( tmax=ee.Number(tmax), tmin=ee.Number(tmin), ea=ee.Number(ea), - fcd=ee.Number(fcd) + fcd=ee.Number(fcd), )) assert float(output) == pytest.approx(expected) @@ -657,7 +657,7 @@ def test_rnl_daily_image(tmin=d_args['tmin'], tmax=d_args['tmax'], ea=d_args['ea'], fcd=d_args['fcd'], expected=d_args['rnl']): output = utils.constant_image_value(calcs._rnl_daily( tmax=ee.Image.constant(tmax), tmin=ee.Number(tmin), - ea=ee.Number(ea), fcd=ee.Number(fcd) + ea=ee.Number(ea), fcd=ee.Number(fcd), )) assert float(output['constant']) == pytest.approx(expected) diff --git a/openet/refetgee/tests/test_daily_refet_output.py b/openet/refetgee/tests/test_daily_refet_output.py index ab24113..a94490d 100644 --- a/openet/refetgee/tests/test_daily_refet_output.py +++ b/openet/refetgee/tests/test_daily_refet_output.py @@ -27,7 +27,7 @@ class DailyData(): csv_df = pd.read_csv(csv_path, engine='python', na_values='NO RECORD') csv_df.rename( columns={'MN': 'TMIN', 'MX': 'TMAX', 'YM': 'TDEW', 'UA': 'WIND', 'SR': 'RS'}, - inplace=True + inplace=True, ) csv_df['DATE'] = csv_df[['YEAR', 'MONTH', 'DAY']].apply( lambda x: dt.datetime(*x).strftime('%Y-%m-%d'), axis=1 @@ -65,11 +65,14 @@ class DailyData(): # Read in the OUT file using pandas (skip header and units) out_df = pd.read_table( out_path, sep='\\s+', index_col=False, - skiprows=list(range(out_start)) + [out_start + 1]) + skiprows=list(range(out_start)) + [out_start + 1], + ) out_df.rename( - columns={'Yr': 'YEAR', 'Mo': 'MONTH', 'Day': 'DAY', 'Tmax': 'TMAX', - 'Tmin': 'TMIN', 'Wind': 'WIND', 'Rs': 'RS', 'DewP': 'TDEW'}, - inplace=True + columns={ + 'Yr': 'YEAR', 'Mo': 'MONTH', 'Day': 'DAY', 'Tmax': 'TMAX', 'Tmin': 'TMIN', + 'Wind': 'WIND', 'Rs': 'RS', 'DewP': 'TDEW', + }, + inplace=True, ) out_df['DATE'] = out_df[['YEAR', 'MONTH', 'DAY']].apply( lambda x: dt.datetime(*x).strftime('%Y-%m-%d'), axis=1 diff --git a/openet/refetgee/tests/test_hourly_refet_output.py b/openet/refetgee/tests/test_hourly_refet_output.py index 4e1c782..4c323cd 100644 --- a/openet/refetgee/tests/test_hourly_refet_output.py +++ b/openet/refetgee/tests/test_hourly_refet_output.py @@ -58,7 +58,8 @@ class HourlyData(): # in2_df.columns = [' '.join(col).replace('-', '').strip() for col in in2_df.columns.values] # in2_df.rename( # columns={'Year': 'YEAR', 'Mo': 'MONTH', 'Da': 'DAY', 'DoY': 'DOY', 'HrMn': 'HOUR'}, - # inplace=True) + # inplace=True, + # ) # in2_df['HOUR'] = (in2_df['HOUR'] / 100).astype(int) # # AgriMet times are local with DST (this will drop one hour) # in2_df['DATETIME'] = in2_df[['YEAR', 'MONTH', 'DAY', 'HOUR']].apply( @@ -76,12 +77,14 @@ class HourlyData(): # Read in the OUT file using pandas (skip header and units) out_df = pd.read_table( out_path, sep='\\s+', index_col=False, - skiprows=list(range(out_start)) + [out_start + 1]) + skiprows=list(range(out_start)) + [out_start + 1], + ) out_df.rename( - columns={'Yr': 'YEAR', 'Mo': 'MONTH', 'Day': 'DAY', 'HrMn': 'HOUR', - 'Tmax': 'TMAX', 'Tmin': 'TMIN', 'DewP': 'TDEW', - 'Wind': 'WIND', 'Rs': 'RS'}, - inplace=True + columns={ + 'Yr': 'YEAR', 'Mo': 'MONTH', 'Day': 'DAY', 'HrMn': 'HOUR', + 'Tmax': 'TMAX', 'Tmin': 'TMIN', 'DewP': 'TDEW', 'Wind': 'WIND', 'Rs': 'RS', + }, + inplace=True, ) out_df['HOUR'] = (out_df['HOUR'] / 100).astype(int) # AgriMet times are local with DST (this will drop one hour) @@ -140,7 +143,7 @@ class HourlyData(): 'lon': lon, # DEADBEEF # 'lat': lat * math.pi / 180, - # 'lon': lon * math.pi / 180 + # 'lon': lon * math.pi / 180, }) values.append(date_values) ids.append('{}-{}'.format(test_date, surface)) diff --git a/openet/refetgee/units.py b/openet/refetgee/units.py index 657e33b..a986d7a 100644 --- a/openet/refetgee/units.py +++ b/openet/refetgee/units.py @@ -4,11 +4,14 @@ def _deg2rad(deg): return deg * math.pi / 180.0 + def _rad2deg(rad): return rad * 180.0 / math.pi + def _c2f(c): return c * (9.0 / 5) + 32 + def _f2c(f): return (f - 32) * (5.0 / 9) diff --git a/pyproject.toml b/pyproject.toml index bb88fea..2ffe822 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "openet-refet-gee" -version = "0.6.3" +version = "0.6.4" authors = [ { name = "Charles Morton", email = "charles.morton@dri.edu" }, ]