package de.hft.stuttgart.water; import java.io.BufferedWriter; import java.io.IOException; import java.nio.charset.StandardCharsets; import java.nio.file.Files; import java.nio.file.Path; import java.nio.file.Paths; import java.time.Duration; import java.time.LocalDateTime; import java.time.format.DateTimeFormatter; import java.util.ArrayList; import java.util.List; import java.util.Locale; import java.util.stream.Collectors; import com.opencsv.CSVReader; import com.opencsv.exceptions.CsvException; public class TreeWaterDemand { /** * Rs/Rso für Nacht = 0.5 * * @param args * @throws IOException * @throws CsvException */ public static void main(String[] args) throws IOException, CsvException { List csvValues; Path csvPath = Paths.get("Wuerzburg-hour-Meteonorm8.csv"); try (CSVReader csvReader = new CSVReader(Files.newBufferedReader(csvPath))) { csvValues = csvReader.readAll(); } int tempIndex = -1; int irrIndex = -1; int humidityIndex = -1; int windSpeedIndex = -1; String[] header = csvValues.get(1); for (int i = 0; i < header.length; i++) { if (header[i].startsWith("GHI (W/m^2)")) { irrIndex = i; } else if (header[i].startsWith("Dry-bulb (C)")) { tempIndex = i; } else if (header[i].startsWith("RHum (%)")) { humidityIndex = i; } else if (header[i].startsWith("Wspd (m/s)")) { windSpeedIndex = i; } } double longitude = 9.933333; // germany = 345 double lz = 345; // Würzburg = double lm = 360 - longitude; double height = 177; double svf = 1; double advection = 1.6; double tr = 0.9; LocalDateTime startOfYear = LocalDateTime.of(2005, 1, 1, 0, 0); LocalDateTime start = LocalDateTime.of(2005, 1, 1, 0, 0); LocalDateTime end = LocalDateTime.of(2006, 1, 1, 0, 0); // LocalDateTime start = LocalDateTime.of(2005, 4, 1, 0, 0); // LocalDateTime end = LocalDateTime.of(2005, 10, 1, 0, 0); List et0us = new ArrayList<>(); List et0s = new ArrayList<>(); List demands = new ArrayList<>(); List rain = new ArrayList<>(); LocalDateTime counter = start; System.out.println("Date,ETIa (L/m²/d),ET0 (L/m²/d)"); int startHours = (int) Duration.between(startOfYear, start).toHours(); int hours = (int) Duration.between(start, end).toHours(); DateTimeFormatter f = DateTimeFormatter.ofPattern("dd/MM/yyyy HH:mm"); Locale.setDefault(new Locale("en", "US")); try (BufferedWriter bw = Files.newBufferedWriter(Paths.get("all_values.csv"), StandardCharsets.UTF_8)) { bw.write("DateTime;Temperature;Humidity;GHI;WindSpeed;SaturationVaporPressure;ActualVaporPressure;ET0;ET0u\n"); bw.write("dd/mm/yyyy HH:MM;[°C];[%];[W/m²];[m/s];[?];[?];[l/h];[l/h]\n"); for (int i = startHours; i < startHours + hours; i++) { bw.write(String.format("%s;", counter.format(f))); String[] row = csvValues.get(i + 2); double temp = Double.parseDouble(row[tempIndex]); double humidity = Double.parseDouble(row[humidityIndex]); double incomingRadiation = Double.parseDouble(row[irrIndex]); bw.write(String.format("%.1f", temp)); bw.write(String.format(";%.0f", humidity)); bw.write(String.format(";%.0f", incomingRadiation)); // convert to MJ/h - 1 W = 0.0036 MJ/h incomingRadiation *= 0.0036; double windSpeed = Double.parseDouble(row[windSpeedIndex]); double saturationVaporPressure = e0(temp); double actualVaporPressure = ea(saturationVaporPressure, humidity); int dayOfYear = counter.getDayOfYear(); double hourOfDay = counter.getHour() + 0.5; double et0 = et0(temp, actualVaporPressure, incomingRadiation, height, windSpeed, humidity, dayOfYear, hourOfDay, lz, lm); double et0u = et0u(svf, advection, et0); bw.write(String.format(";%.1f", windSpeed)); bw.write(String.format(";%.3f", saturationVaporPressure)); bw.write(String.format(";%.3f", actualVaporPressure)); bw.write(String.format(";%.3f", et0)); bw.write(String.format(";%.3f", et0u)); et0s.add(et0); // System.out.println("et0: " + et0); // System.out.println("etia: " + etia); et0us.add(et0u); counter = counter.plusHours(1); bw.write("\n"); } } // System.out.println(rain.stream().collect(Collectors.summarizingDouble(Double::doubleValue))); Double et0Sum = et0s.stream().collect(Collectors.summingDouble(Double::doubleValue)); System.out.println("Summe aller ET0: " + et0Sum); Double et0uSum = et0us.stream().collect(Collectors.summingDouble(Double::doubleValue)); System.out.println("Summe aller ET0u: " + et0uSum); double et0uFactor = (0.865 * Math.log10(1 / et0uSum) + 3.36); for (double et0u : et0us) { double etia = 0; if (et0u > 0) { etia = etia(tr, et0u, et0uFactor); } demands.add(etia); } System.out .println("Summe aller ETIa: " + demands.stream().collect(Collectors.summingDouble(Double::doubleValue))); // printHourly(start, et0s, demands); // printDaily(start, et0s, demands); // test(); } private static void printHourly(LocalDateTime start, List et0s, List demands) { LocalDateTime date = start; for (int i = 0; i < demands.size(); i++) { System.out.println(date + "," + demands.get(i) + "," + et0s.get(i)); date = date.plusHours(1); } } private static void printDaily(LocalDateTime start, List et0s, List demands) { LocalDateTime date = start; for (int i = 0; i < demands.size(); i = i + 24) { double demandSum = 0; double et0Sum = 0; for (int j = i; j < i + 24; j++) { demandSum += demands.get(j); et0Sum += et0s.get(j); } System.out.println(date + "," + demandSum + "," + et0Sum); date = date.plusHours(24); } } private static void test() { double lz = 15; double lm = 16.25; double height = 8; double temp = 38; double humidity = 52; double windSpeed = 3.3; double incomingRadiation = 2.45; double saturationVaporPressure = e0(temp); double actualVaporPressure = ea(saturationVaporPressure, humidity); int dayOfYear = 274; double hourOfDay = 14 + 0.5; double et0 = et0(temp, actualVaporPressure, incomingRadiation, height, windSpeed, humidity, dayOfYear, hourOfDay, lz, lm); System.out.println(et0); } private static double etia(double tr, double et0u, double et0uFactor) { // return 1 * 1.6 * et0u; return tr * 1.25 * et0u * et0uFactor; // return tr * (1.61 * Math.log10(800) - 3.39) * et0u * (0.865 * Math.log10(1 / et0u) + 3.36); } private static double et0u(double svf, double advection, double et0) { return svf * advection * et0; } private static double et0(double temperature, double actualVaporPressure, double shortWaveRadiation, double height, double windSpeed, double humidity, int dayOfYear, double hourOfDay, double lz, double lm) { double nominalRadiation = rn(temperature, actualVaporPressure, shortWaveRadiation, height, dayOfYear, hourOfDay, lz, lm); double gamma = gamma(height); double saturationVaporPressure = e0(temperature); double delta = delta(temperature, saturationVaporPressure); double g = g(nominalRadiation, shortWaveRadiation > 0); double n = 0.408 * delta * (nominalRadiation - g) + gamma * (37 / (temperature + 273)) * windSpeed * (saturationVaporPressure - actualVaporPressure); return n / (delta + gamma * (1 + 0.34 * windSpeed)); } private static double ea(double saturationVaporPressure, double humidity) { return saturationVaporPressure * humidity / 100; } private static double e0(double temperature) { return 0.6108 * Math.exp((17.27 * temperature) / (temperature + 237.3)); } private static double gamma(double height) { double pressure = 101.3 * Math.pow((293 - 0.0065 * height) / 293, 5.26); return 0.665 * 0.001 * pressure; } private static double g(double nominalRadiation, boolean daylight) { if (daylight) { return 0.1 * nominalRadiation; } else { return 0.5 * nominalRadiation; } } private static double rn(double temperature, double actualVaporPressure, double shortWaveRadiation, double height, int dayOfYear, double hourOfDay, double lz, double lm) { // 1 - 0.23: albedo effect of grass return (1 - 0.23) * shortWaveRadiation - rnl(temperature, actualVaporPressure, shortWaveRadiation, height, dayOfYear, hourOfDay, lz, lm); } private static double rnl(double temperature, double actualVaporPressure, double shortWaveRadiation, double height, int dayOfYear, double hourOfDay, double lz, double lm) { // divide by 24 for hourly double sigma = 4.903E-9 / 24; double rs = shortWaveRadiation; double dr = 1 + 0.033 * Math.cos(2 * Math.PI * dayOfYear / 365d); double b = (2 * Math.PI * (dayOfYear - 81)) / 364; double sc = 0.1645 * Math.sin(2 * b) - 0.1255 * Math.cos(b) - 0.025 * Math.sin(b); double omega = Math.PI / 12 * ((hourOfDay + 0.06667 * (lz - lm) + sc) - 12); double phi = lm * Math.PI / 180; double delta = 0.409 * Math.sin(2 * Math.PI / 365 * dayOfYear - 1.39); double omegas = Math.acos(-Math.tan(phi) * Math.tan(delta)); double ra = 0; if (!(omega < -omegas || omega > omegas)) { double omega1 = omega - Math.PI * 1 / 24; double omega2 = omega + Math.PI * 1 / 24; ra = 12 * 60 / Math.PI * 0.0820 * dr * ((omega2 - omega1) * Math.sin(phi) * Math.sin(delta) + Math.cos(phi) * Math.cos(delta) * (Math.sin(omega2) - Math.sin(omega1))); } double rso = (0.75 + 2 * 10E-5 * height) * ra; double rsRsoQuotient = 0.5; if (rso != 0) { rsRsoQuotient = rs / rso; } if (rsRsoQuotient > 1) { rsRsoQuotient = 1; } return sigma * Math.pow(temperature + 273.16, 4) * (0.34 - 0.14 * Math.sqrt(actualVaporPressure)) * (1.35 * rsRsoQuotient - 0.35); } private static double delta(double temperatur, double e0) { return (4098 * e0) / Math.pow(temperatur + 237.3, 2); } }