Tentative solution to handle daily radiation properly - also fixes for same matter in ET

Author: moovida

hmachine/src/main/java/org/hortonmachine/hmachine/geoframe/ermworkflow/ErmPrestleyEt.java

@@ -9,6 +9,7 @@
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.BasinDataSchema.BasinDataField;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.EnvironmentalVariableType;
+import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.TimeResolution;
 import org.hortonmachine.hmachine.geoframe.utils.IWaterBudgetSimulationRunner;
 import org.hortonmachine.hmachine.geoframe.utils.PrestleyETAtCentroid;
 
@@ -64,12 +65,21 @@ public class ErmPrestleyEt extends HMModel {
 	@In
 	public String pEndTimestamp;
 
+	@Description("The expected time resolution of the data. Daily and hourly (default) is supported.")
+	@In
+	public TimeResolution pTimeResolution = TimeResolution.HOURLY;
+
 	@Description("If true, existing output files are overwritten.")
 	@In
 	public boolean doOverwrite = false;
 
 	@Execute
 	public void process() throws Exception {
+		if (pTimeResolution == TimeResolution.MONTHLY || pTimeResolution == TimeResolution.YEARLY) {
+			throw new UnsupportedOperationException(
+					"ErmPrestleyEt only supports HOURLY and DAILY resolutions, got " + pTimeResolution);
+		}
+
 		try (ASpatialDb db = EDb.GEOPACKAGE.getSpatialDb()) {
 			db.open(inGpkg);
 
@@ -99,7 +109,7 @@ public void process() throws Exception {
 
 			var ptEt = new PrestleyETAtCentroid();
 			ptEt.inGeoframeDB = db;
-			ptEt.isHourly = true;
+			ptEt.isHourly = pTimeResolution == TimeResolution.HOURLY;
 			ptEt.pAlpha = 1.26;
 			ptEt.inTempReader = temperatureReader;
 			ptEt.inNetReader = netReader;

hmachine/src/main/java/org/hortonmachine/hmachine/geoframe/ermworkflow/ErmRadiation.java

@@ -42,6 +42,16 @@
  * database; humidity and the atmospheric clearness index are left at their
  * clear-sky defaults, since this launcher assumes clear-sky conditions with
  * only temperature known.
+ *
+ * <p>
+ * Note that the output is always a flux, in W/m2: 
+ * <ul>
+ * <li>for {@link TimeResolution#HOURLY} it is the instantaneous net
+ * radiation for that hour</li>
+ * <li>for {@link TimeResolution#DAILY} - the
+ * avg of 24 hourly samples taken at that day's actual sun positions
+ * </li>
+ * </ul>
  */
 @Description("Radiation calculator.")
 @Author(name = "Daniele Andreis", contact = "")
@@ -69,10 +79,15 @@ public class ErmRadiation extends HMModel {
 	@In
 	public String pEndTimestamp;
 
-	@Description("Time resolution. Defaults to hourly")
+	@Description("The expected time resolution of the data. Daily and hourly (default) is supported.")
 	@In
 	public TimeResolution pTimeResolution = TimeResolution.HOURLY;
 
+	@Description("Number of sun-position samples used to average net radiation over a day if "
+			+ "pTimeResolution is DAILY. 24 (one per hour) is the most accurate but also slowest.")
+	@In
+	public int pDailySubSamples = 24;
+
 	@Description("If true, existing output files are overwritten.")
 	@In
 	public boolean doOverwrite = false;
@@ -83,6 +98,11 @@ public class ErmRadiation extends HMModel {
 
 	@Execute
 	public void process() throws Exception {
+		if (pTimeResolution == TimeResolution.MONTHLY || pTimeResolution == TimeResolution.YEARLY) {
+			throw new UnsupportedOperationException(
+					"ErmRadiation only supports HOURLY and DAILY resolutions, got " + pTimeResolution);
+		}
+
 		Paths p = new Paths(inDtm, doOverwrite);
 
 		try (ASpatialDb db = EDb.GEOPACKAGE.getSpatialDb()) {
@@ -116,7 +136,8 @@ public void process() throws Exception {
 			radiation.dem =  dtm; // TODO Daniele, why where you using the pit here?
 			radiation.inSkyview = skyview;
 			radiation.lwrvModeel = "6";
-			radiation.doHourly = pTimeResolution == TimeResolution.HOURLY;
+			radiation.pTimeResolution = pTimeResolution;
+			radiation.pDailySubSamples = pDailySubSamples;
 			radiation.init();
 			radiation.process();
 

hmachine/src/main/java/org/hortonmachine/hmachine/geoframe/io/database/TableUtils.java

@@ -6,9 +6,6 @@
 import org.hortonmachine.dbs.compat.ASpatialDb;
 import org.hortonmachine.dbs.compat.objects.QueryResult;
 import org.hortonmachine.gears.libs.modules.HMConstants;
-import org.hortonmachine.hmachine.geoframe.io.database.tables.GeoFrameGeoTable;
-import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.StationSchema.Station;
-import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.StationSchema.StationType;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.EnvironmentalVariable;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.EnvironmentalVariableType;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.TimeResolution;
@@ -26,27 +23,17 @@ public final static List<EnvironmentalVariable> getFixedEnviramentalVariable(Tim
 
 		String mmFlux = "mm/";
 		String temperatureUnit = "°C";
-		String radiationUnit = null;
+		
+		// Radiation is always a flux and in W/m² and is always related to the size of the timestep.
+		String radiationUnit = "W/m²";
 		String dischargeUnit = "m³/s";
 
 		if (resolution != null) {
 			switch (resolution) {
-			case HOURLY -> {
-				mmFlux = mmFlux + "h";
-				radiationUnit = "W/m²";
-			}
-			case DAILY -> {
-				mmFlux = mmFlux + "day";
-				radiationUnit = "MJ/m²/day";
-			}
-			case MONTHLY -> {
-				mmFlux = mmFlux + "month";
-				radiationUnit = "MJ/m²/month";
-			}
-			case YEARLY -> {
-				mmFlux = mmFlux + "year";
-				radiationUnit = "MJ/m²/year";
-			}
+			case HOURLY -> mmFlux = mmFlux + "h";
+			case DAILY -> mmFlux = mmFlux + "day";
+			case MONTHLY -> mmFlux = mmFlux + "month";
+			case YEARLY -> mmFlux = mmFlux + "year";
 			}
 		} else {
 			mmFlux = null;

hmachine/src/main/java/org/hortonmachine/hmachine/geoframe/utils/RadiationAtCentroid.java

@@ -27,6 +27,7 @@
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.BasinPolygonSchema;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.BasinPolygonSchema.BasinMultiPolygonField;
 import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.EnvironmentalVariableType;
+import org.hortonmachine.hmachine.geoframe.io.database.tables.implementation.VarSchema.TimeResolution;
 import org.hortonmachine.hmachine.geoframe.utils.radiation.NetRadiationPointCase;
 import org.hortonmachine.hmachine.geoframe.utils.radiation.LwrbPointCase.Lwrb;
 import org.hortonmachine.hmachine.geoframe.utils.radiation.swrbPointCase.ShortwaveRadiationBalancePointCase;
@@ -80,10 +81,15 @@ public class RadiationAtCentroid extends HMModel {
 	@In
 	public double bCloud = 1;
 
-	@Description("Toggle to compute hourly radiation values, if false daily values are computed")
+	@Description("The expected time resolution of the data. Daily and hourly (default) is supported.")
 	@In
-	public boolean doHourly = true;
-	
+	public TimeResolution pTimeResolution = TimeResolution.HOURLY;
+
+	@Description("Number of sun-position samples used to average net radiation over a day if "
+			+ "pTimeResolution is DAILY. 24 (one per hour) is the most accurate but also slowest.")
+	@In
+	public int pDailySubSamples = 24;
+
 	public double alpha = 0.26;
 
 	@Description("Ozone layer thickness in cm")
@@ -129,11 +135,21 @@ public class RadiationAtCentroid extends HMModel {
 	private HashMap<Integer, double[]> inNan;
 	private SimpleFeatureCollection inBasinsFC;
 
+	private static final long MILLIS_PER_HOUR = 3_600_000L;
+	private static final long MILLIS_PER_DAY = 24 * MILLIS_PER_HOUR;
+
 	@Initialize
 	public void init() throws Exception {
 		if (lwrvModeel == null || lwrvModeel.isEmpty()) {
 			throw new IllegalArgumentException();
 		}
+		if (pTimeResolution == TimeResolution.MONTHLY || pTimeResolution == TimeResolution.YEARLY) {
+			throw new UnsupportedOperationException(
+					"ErmRadiation only supports HOURLY and DAILY resolutions, got " + pTimeResolution);
+		}
+		if (pDailySubSamples < 1) {
+			throw new IllegalArgumentException("pDailySubSamples must be >= 1, got " + pDailySubSamples);
+		}
 		checkNull(inGeoframeDb);
 
 		if (!(inGeoframeDb.hasTable(GeoFrameGeoTable.BASIN.tableName())
@@ -201,7 +217,7 @@ private Lwrb createLwrb() {
 	 */
 	private ShortwaveRadiationBalancePointCase createSwrb() {
 		ShortwaveRadiationBalancePointCase s = new ShortwaveRadiationBalancePointCase();
-		s.doHourly = doHourly;
+		s.doHourly = pTimeResolution == TimeResolution.HOURLY;
 		s.fStationsid = BasinMultiPolygonField.ID.columnName();
 		s.inStationsFC = inBasinsFC;
 		s.inDem = dem;
@@ -289,23 +305,8 @@ public void process() throws Exception {
 
 	private void processTimestep(HashMap<Integer, double[]> temperature, HashMap<Integer, double[]> humidity,
 			HashMap<Integer, double[]> clearSky, long t) throws Exception {
-		lwrb.inAirTemperatureValuesHM = temperature;
-		lwrb.inSoilTempratureValuesHM = temperature;
-		lwrb.inHumidityValuesHM = humidity;
-		lwrb.inClearnessIndexValuesHM = clearSky;
-		swrb.inHumidityValues = humidity;
-		swrb.inTemperatureValues = temperature;
-		swrb.tCurrentDateString = GeoframeEnvDatabaseIterator.ts2str(t);
-		
-		lwrb.process();
-		swrb.process();
-		
-		nrpc.inShortwaveValues = swrb.outHMtotal;
-		nrpc.inDownwellingValues = lwrb.outHMlongwaveDownwellingHM;
-		nrpc.inUpwellingValues = lwrb.outHMlongwaveUpwellingHM;
-		nrpc.process();
-		
-		HashMap<Integer, double[]> out = nrpc.outHMnetRad;
+		HashMap<Integer, double[]> out = computeNetRadiation(lwrb, swrb, nrpc, temperature, humidity, clearSky, t);
+
 		String insertSql = GeoFrameSimpleTable.BASINDATA.getSchema().buildInsertAll();
 		inGeoframeDb.execOnConnection(conn -> {
 			boolean autoCommit = conn.getAutoCommit();
@@ -350,9 +351,7 @@ private static class TimestepResult {
 	 * <p>
 	 * The result map is copied out of {@code nrpc.outHMnetRad} because that
 	 * field is reused and overwritten in place on every {@code process()}
-	 * call on the same instance - if a thread goes on to compute another
-	 * timestep before this result is consumed, an uncopied reference would
-	 * end up pointing at that later timestep's data instead.
+	 * call on the same instance.
 	 */
 	private TimestepResult computeTimestepParallel(int idx, int[] ids, ThreadLocal<Lwrb> tlLwrb,
 			ThreadLocal<ShortwaveRadiationBalancePointCase> tlSwrb, ThreadLocal<NetRadiationPointCase> tlNrpc)
@@ -376,23 +375,95 @@ private TimestepResult computeTimestepParallel(int idx, int[] ids, ThreadLocal<L
 		ShortwaveRadiationBalancePointCase swrbLocal = tlSwrb.get();
 		NetRadiationPointCase nrpcLocal = tlNrpc.get();
 
+		HashMap<Integer, double[]> out = computeNetRadiation(lwrbLocal, swrbLocal, nrpcLocal, temperature, humidity,
+				clearSky, timestamp);
+		return new TimestepResult(timestamp, out);
+	}
+
+	/**
+	 * Computes net radiation (W/m2) for one timestep, dispatching to
+	 * either a single hourly sample or a daily average of 24 hourly samples
+	 * depending on the time resolution.
+	 *
+	 * <p>
+	 * The returned map is always a fresh copy, independent of
+	 * {@code nrpc.outHMnetRad}, which is reused and overwritten in place on
+	 * every {@code process()} call - without a copy, a thread going on to
+	 * compute another timestep before this result is consumed would leave the
+	 * caller holding a reference to that later timestep's data instead.
+	 */
+	private HashMap<Integer, double[]> computeNetRadiation(Lwrb lwrbLocal, ShortwaveRadiationBalancePointCase swrbLocal,
+			NetRadiationPointCase nrpcLocal, HashMap<Integer, double[]> temperature, HashMap<Integer, double[]> humidity,
+			HashMap<Integer, double[]> clearSky, long t) throws Exception {
 		lwrbLocal.inAirTemperatureValuesHM = temperature;
 		lwrbLocal.inSoilTempratureValuesHM = temperature;
 		lwrbLocal.inHumidityValuesHM = humidity;
 		lwrbLocal.inClearnessIndexValuesHM = clearSky;
+		// Longwave depends only on temperature/humidity/clearness - not on time of day -
+		// so one process() call covers the whole timestep, hourly or daily alike.
+		lwrbLocal.process();
+
 		swrbLocal.inHumidityValues = humidity;
 		swrbLocal.inTemperatureValues = temperature;
-		swrbLocal.tCurrentDateString = GeoframeEnvDatabaseIterator.ts2str(timestamp);
 
-		lwrbLocal.process();
-		swrbLocal.process();
+		if (pTimeResolution == TimeResolution.HOURLY) {
+			swrbLocal.doHourly = true;
+			swrbLocal.tCurrentDateString = GeoframeEnvDatabaseIterator.ts2str(t);
+			swrbLocal.process();
 
-		nrpcLocal.inShortwaveValues = swrbLocal.outHMtotal;
-		nrpcLocal.inDownwellingValues = lwrbLocal.outHMlongwaveDownwellingHM;
-		nrpcLocal.inUpwellingValues = lwrbLocal.outHMlongwaveUpwellingHM;
-		nrpcLocal.process();
+			nrpcLocal.inShortwaveValues = swrbLocal.outHMtotal;
+			nrpcLocal.inDownwellingValues = lwrbLocal.outHMlongwaveDownwellingHM;
+			nrpcLocal.inUpwellingValues = lwrbLocal.outHMlongwaveUpwellingHM;
+			nrpcLocal.process();
+			return new HashMap<>(nrpcLocal.outHMnetRad);
+		} else {
+			return averageDailyNetRadiation(lwrbLocal, swrbLocal, nrpcLocal, t);
+		}
+	}
 
-		return new TimestepResult(timestamp, new HashMap<>(nrpcLocal.outHMnetRad));
+	/**
+	 * Averages net radiation over a full day into a daily-mean flux.
+	 *
+	 * <p>
+	 * {@link ShortwaveRadiationBalancePointCase} only ever evaluates the sun's
+	 * position (and the DEM-derived shading/skyview weighting that goes with
+	 * it) at one instant, and that DEM-wide shadow recomputation is what
+	 * dominates runtime. With daily-resolution input there is a single
+	 * temperature/humidity/clearness reading for the whole day, so this
+	 * processes that same daily reading at {@link #pDailySubSamples} sun
+	 * positions evenly spaced across the day and averages the results. At the
+	 * default of 24 this means running hourly resolution for the day and
+	 * averaging the 24 outputs; lower values are cheaper but coarser
+	 * approximations, useful for quick test runs.
+	 */
+	private HashMap<Integer, double[]> averageDailyNetRadiation(Lwrb lwrbLocal,
+			ShortwaveRadiationBalancePointCase swrbLocal, NetRadiationPointCase nrpcLocal, long dayTimestamp)
+			throws Exception {
+		long dayStartMillis = Math.floorDiv(dayTimestamp, MILLIS_PER_DAY) * MILLIS_PER_DAY;
+		long stepMillis = MILLIS_PER_DAY / pDailySubSamples;
+
+		HashMap<Integer, Double> sumWattsPerM2 = new HashMap<>();
+		swrbLocal.doHourly = true; // sample the actual sun position each time, rather than one fixed instant
+		for (int sample = 0; sample < pDailySubSamples; sample++) {
+			long sampleMillis = dayStartMillis + sample * stepMillis;
+			swrbLocal.tCurrentDateString = GeoframeEnvDatabaseIterator.ts2str(sampleMillis);
+			swrbLocal.process();
+
+			nrpcLocal.inShortwaveValues = swrbLocal.outHMtotal;
+			nrpcLocal.inDownwellingValues = lwrbLocal.outHMlongwaveDownwellingHM;
+			nrpcLocal.inUpwellingValues = lwrbLocal.outHMlongwaveUpwellingHM;
+			nrpcLocal.process();
+
+			for (Map.Entry<Integer, double[]> entry : nrpcLocal.outHMnetRad.entrySet()) {
+				sumWattsPerM2.merge(entry.getKey(), entry.getValue()[0], Double::sum);
+			}
+		}
+
+		HashMap<Integer, double[]> dailyMean = new HashMap<>();
+		for (Map.Entry<Integer, Double> entry : sumWattsPerM2.entrySet()) {
+			dailyMean.put(entry.getKey(), new double[] { entry.getValue() / pDailySubSamples });
+		}
+		return dailyMean;
 	}
 
 	/**

hmachine/src/main/java/org/hortonmachine/hmachine/modules/hydrogeomorphology/etp/OmsPresteyTaylorEtpModel.java

@@ -167,7 +167,7 @@ public void process() throws Exception {
 			}
 
 			double etp = getET(pGmorn, pGnight, pAlpha, inputNetRadiation,
-					doHourly ? defaultHourlyNetradiation : defaultHourlyNetradiation, t, defaultTemp, pressure,
+					doHourly ? defaultHourlyNetradiation : defaultDailyNetradiation, t, defaultTemp, pressure,
 					defaultPressure, doHourly, tCurrent);
 			outPTEtp.put(basinId, new double[] { etp });
 		}