Prototype

.classpath

@@ -0,0 +1,9 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<classpath>
+	<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER/org.eclipse.jdt.internal.debug.ui.launcher.StandardVMType/JavaSE-1.8"/>
+	<classpathentry kind="src" path="src"/>
+	<classpathentry kind="src" path="tests"/>
+	<classpathentry kind="con" path="org.eclipse.jdt.junit.JUNIT_CONTAINER/5"/>
+	<classpathentry kind="lib" path="libs/tomlj-1.1.0-all.jar"/>
+	<classpathentry kind="output" path="bin"/>
+</classpath>

.github/workflows/tests.yml

@@ -0,0 +1,20 @@
+name: Run Tests
+
+on:
+  push:
+    branches:
+      - main
+
+jobs:
+  build-and-test:
+    name: Build and test
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v3
+    - uses: actions/setup-java@v3
+      with:
+        distribution: 'zulu'
+        java-version: '8'
+        java-package: jdk
+        architecture: x64
+    - run: ci/test.sh

.gitignore

@@ -0,0 +1,3 @@
+*.bak
+*.class
+.settings

.project

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<projectDescription>
+	<name>SimRasterise</name>
+	<comment></comment>
+	<projects>
+	</projects>
+	<buildSpec>
+		<buildCommand>
+			<name>org.eclipse.jdt.core.javabuilder</name>
+			<arguments>
+			</arguments>
+		</buildCommand>
+	</buildSpec>
+	<natures>
+		<nature>org.eclipse.jdt.core.javanature</nature>
+	</natures>
+</projectDescription>

ci/test.bat

@@ -0,0 +1,3 @@
+md build
+javac -classpath .;libs/junit-platform-console-standalone-1.9.2.jar;tomlj-1.1.0-all.jar -d build src/com/mrc/rasterise/*.java src/com/mrc/rasterise/population/*.java src/com/mrc/rasterise/shapefile/*.java
+java -jar libs/junit-platform-console-standalone-1.9.2.jar --class-path build --scan-class-path

ci/test.sh

@@ -0,0 +1,3 @@
+mkdir -p build
+javac -classpath .:libs/junit-platform-console-standalone-1.9.2.jar -d build src/rasterise/*.java
+java -jar libs/junit-platform-console-standalone-1.9.2.jar --class-path build --scan-class-path

src/com/mrc/rasterise/Config.java

@@ -0,0 +1,41 @@
+package com.mrc.rasterise;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+import java.nio.file.Paths;
+
+import org.tomlj.Toml;
+import org.tomlj.TomlParseResult;
+
+public class Config  {
+
+  TomlParseResult toml;
+
+  void validate() throws Exception {
+    if (!toml.contains("population")) {
+      throw new Exception("population not found in config");
+    }
+
+    if (!new File(toml.getString("population")).exists()) {
+      throw new FileNotFoundException("Population file not found : "+toml.getString("population"));
+    }
+  }
+
+  public Config(String file) throws Exception {
+    if (!new File(file).exists()) {
+      throw new FileNotFoundException("Config file not found: "+file);
+    }
+
+    toml = Toml.parse(Paths.get(file));
+
+    int errorCount = toml.errors().size();
+    if (errorCount > 0) {
+      StringBuilder sb = new StringBuilder ();
+      toml.errors().forEach(error -> sb.append(sb+"\n"));
+      throw new Exception("Errors in config file: \n"+ sb.toString());
+    }
+
+    validate();
+  }
+}
+

src/com/mrc/rasterise/Rasterise.java

@@ -0,0 +1,26 @@
+package com.mrc.rasterise;
+
+import java.io.InputStream;
+
+public class Rasterise {
+
+  public int check_args(String[] args) throws Exception {
+    if (args.length != 1) {
+      throw new Exception("Requires one argument: config.toml");
+    }
+    return 0;
+  }
+
+  public int start(InputStream input, String[] args) throws Exception {
+    check_args(args);
+    Config config = new Config(args[0]);
+    System.out.println(config.toml.get("population"));
+    return 0;
+  }
+
+  public static void main(String[] args) throws Exception {
+    Rasterise r = new Rasterise();
+    r.start(System.in, args);
+  }
+
+}

src/com/mrc/rasterise/era5_land/AdminUnitStats.java

@@ -0,0 +1,76 @@
+package com.mrc.rasterise.era5_land;
+
+public class AdminUnitStats {
+  double mint_lo, mint_avg, mint_hi, meant_lo, meant_avg, meant_hi, maxt_lo, maxt_avg, maxt_hi;
+  double prec_lo, prec_avg, prec_hi;
+  double minr_lo, minr_avg, minr_hi, meanr_lo, meanr_avg, meanr_hi, maxr_lo, maxr_avg, maxr_hi;
+
+  int count;
+  void init() {
+    count = 0;
+    mint_lo = Double.POSITIVE_INFINITY;
+    mint_avg = 0;
+    mint_hi = Double.NEGATIVE_INFINITY;
+    meant_lo = Double.POSITIVE_INFINITY;
+    meant_avg = 0;
+    meant_hi = Double.NEGATIVE_INFINITY;
+    maxt_lo = Double.POSITIVE_INFINITY;
+    maxt_avg = 0;
+    maxt_hi = Double.NEGATIVE_INFINITY;
+    minr_lo = Double.POSITIVE_INFINITY;
+    minr_avg = 0;
+    minr_hi = Double.NEGATIVE_INFINITY;
+    meanr_lo = Double.POSITIVE_INFINITY;
+    meanr_avg = 0;
+    meanr_hi = Double.NEGATIVE_INFINITY;
+    maxr_lo = Double.POSITIVE_INFINITY;
+    maxr_avg = 0;
+    maxr_hi = Double.NEGATIVE_INFINITY;
+    prec_lo = Double.POSITIVE_INFINITY;
+    prec_avg = 0;
+    prec_hi = Double.NEGATIVE_INFINITY;
+
+  }
+
+  AdminUnitStats() {
+    init();
+  }
+
+  void update(OneDay D, int x, int y, int pop) {
+    count+=pop;
+    mint_lo = Math.min(mint_lo, D.tm_min[x][y]);
+    mint_avg += (D.tm_min[x][y] * pop);
+    mint_hi = Math.max(mint_hi, D.tm_min[x][y]);
+    meant_lo = Math.min(meant_lo, D.tm_mean[x][y]);
+    meant_avg += (D.tm_mean[x][y] * pop);
+    meant_hi = Math.max(meant_hi, D.tm_mean[x][y]);
+    maxt_lo = Math.min(maxt_lo, D.tm_max[x][y]);
+    maxt_avg += (D.tm_max[x][y] * pop);
+    maxt_hi = Math.max(maxt_hi, D.tm_max[x][y]);
+
+    minr_lo = Math.min(minr_lo, D.rh_min[x][y]);
+    minr_avg += (D.rh_min[x][y] * pop);
+    minr_hi = Math.max(minr_hi, D.rh_min[x][y]);
+    meanr_lo = Math.min(meanr_lo, D.rh_mean[x][y]);
+    meanr_avg += (D.rh_mean[x][y] * pop);
+    meanr_hi = Math.max(meanr_hi, D.rh_mean[x][y]);
+    maxr_lo = Math.min(maxr_lo, D.rh_max[x][y]);
+    maxr_avg += (D.rh_max[x][y] * pop);
+    maxr_hi = Math.max(maxr_hi, D.rh_max[x][y]);
+
+    prec_lo = Math.min(prec_lo, D.precip[x][y]);
+    prec_avg += (D.precip[x][y] * pop);
+    prec_hi = Math.max(prec_hi, D.precip[x][y]);
+
+  }
+
+  void finalise() {
+    mint_avg /= (double) count;
+    meant_avg /= (double) count;
+    maxt_avg /= (double) count;
+    prec_avg /= (double) count;
+    minr_avg /= (double) count;
+    meanr_avg /= (double) count;
+    maxr_avg /= (double) count;
+  }
+}

src/com/mrc/rasterise/era5_land/Coords.java

@@ -0,0 +1,59 @@
+package com.mrc.rasterise.era5_land;
+
+public class Coords {
+
+  public static int lon_to_ls(float lon) {
+    return (int) Math.floor((lon + 180.0) * 120.0);
+  }
+
+  public static int lat_to_ls(float lat) {
+    return (int) Math.floor((90 - lat) * 120.0);
+
+  }
+
+  public static int[] ls_to_era5_x() {
+    // A lookup table to get around some awkward fencepost maths.
+    // Landscan x=0 is the grid-square with longitude -180.0 ... -179.9917
+    // ERA5 x=1800 is the point sample at -180,0.
+
+    // Landscan cells with x = 0, 1, 2, 3, 4 and 5 are all closest to ERA5 where x = 1800.
+    // Landscan cells with x = 6, 7, 8 ,9, 10, 11 are closest to ERA5 x=1801 - as are Landscan cells 12, 13, 14, 15, 16, 17.
+    // Then landscan 18 is closest to 1802....
+
+    // This continues all the way to ERA5 x=3599, which wraps to x=0, and then all the way to x=1799.
+    // Finally, we have the last 6 cells on the extreme east of Landscan's world, which are closest to x=1800 in ERA5 world.
+
+    int[] conv = new int[43200];
+    for (int i = 0; i <= 5; i++) conv[i] = 1800;
+    for (int i = 43199; i >= 43194; i--) conv[i] = 1800;
+    int era5_x = 1801;
+    for (int i = 6; i < 43194; i+=12) {
+      for (int ii = i; ii < i + 12; ii++) {
+        conv[ii] = era5_x;
+      }
+      era5_x = (era5_x + 1) % 3600;
+    }
+    return conv;
+  }
+
+  public static int[] ls_to_era5_y() {
+    // Similarly...ish
+    // Landscan y=0 is the grid-square with lat +89.9917 ... +90
+    // ERA5 y = 0 is at +90 - the top edge.
+    // Landscan y = 0, 1, 2, 3, 4, 5 are closest to ERA5 y = 0
+    // Landscan y = 6, 7, 8 , 9, 10, 11 are closest to ERA5 y = 1, as are ls_y = 12,13,14,15,16,17
+    // At the bottom, Landscan y = 21594..21599 are closest to ERA5 y = 1801, which is on the -90 gridline.
+
+    int[] conv = new int[21600];
+    for (int j = 0; j <= 5; j++) conv[j] = 0;
+    for (int j = 21599; j >= 21594; j--) conv[j] = 1801;
+    int era5_y = 1;
+    for (int j = 6; j < 21594; j+=12) {
+      for (int jj = j; jj < j + 12; jj++) {
+        conv[jj] = era5_y;
+      }
+      era5_y++;
+    }
+    return conv;
+  }
+}

src/com/mrc/rasterise/era5_land/Download.java

@@ -0,0 +1,80 @@
+package com.mrc.rasterise.era5_land;
+
+import java.io.File;
+import java.io.FileWriter;
+import java.io.PrintWriter;
+import java.util.GregorianCalendar;
+
+public class Download {
+
+  private static String leadingZero(int m) {
+    return ((m<10) ? "0" : "") + String.valueOf(m);
+  }
+
+
+  private static void writeDownload(String variable, String year, String month, String day, String time, int proc,
+                                    String file, String workingDir) throws Exception {
+
+    PrintWriter PW = new PrintWriter(new File(workingDir + "fetch"+proc+".py"));
+    PW.println("import cdsapi");
+    PW.println("c = cdsapi.Client()");
+    PW.println("c.retrieve('reanalysis-era5-land', {");
+    PW.println("             'variable': "+variable+",");
+    PW.println("             'year' : '"+year+"',");
+    PW.println("             'month' : '"+month+"',");
+    PW.println("             'day' : '"+day+"',");
+    if (time != null) PW.println("             'time' : "+time+",");
+    PW.println("             'format': 'grib'");
+    PW.println("            }, '"+workingDir+file+"')");
+    PW.close();
+  }
+
+  static Process[] processes = new Process[3];
+
+  private static void pull_data(int min_year, int max_year, String workingDir) throws Exception {
+    GregorianCalendar gc = new GregorianCalendar();
+    gc.setTimeInMillis(System.currentTimeMillis());
+    gc.set(GregorianCalendar.DAY_OF_MONTH,  15);
+
+
+    for (int year = min_year; year <= max_year; year++) {
+      gc.set(GregorianCalendar.YEAR, year);
+      for (int month = 1; month <= 12; month++) {
+        gc.set(GregorianCalendar.DAY_OF_MONTH,  1);
+        gc.set(GregorianCalendar.MONTH, month - 1);
+        int no_days = gc.getActualMaximum(GregorianCalendar.DAY_OF_MONTH);
+        for (int day = 1; day <= no_days; day++) {
+          gc.set(GregorianCalendar.DAY_OF_MONTH, day);
+          for (int type = 0; type <= 1; type++) {
+            int proc = -1;
+            while (proc == -1) {
+              for (int p = 0; p < processes.length; p++) {
+                if ((processes[p] == null) || (!processes[p].isAlive())) {
+                  proc = p;
+                  break;
+                }
+              }
+              if (proc == -1) Thread.sleep(10000);
+            }
+
+            if (type == 0) {
+              writeDownload("'total_precipitation'", String.valueOf(year), leadingZero(month), leadingZero(day), "'00:00'",
+                  proc, "precip_" + year + leadingZero(month)+leadingZero(day) + ".grib", workingDir);
+            } else {
+              writeDownload("['2m_temperature', '2m_dewpoint_temperature']", String.valueOf(year), leadingZero(month), leadingZero(day),
+                  "['00:00', '01:00', '02:00', '03:00', '04:00', '05:00', '06:00', '07:00', '08:00', '09:00', '10:00', '11:00', "+
+                  "'12:00', '13:00', '14:00', '15:00', '16:00', '17:00', '18:00', '19:00', '20:00', '21:00', '22:00', '23:00']", proc,
+                  "temp_" + year + leadingZero(month) + leadingZero(day) + ".grib", workingDir);
+            }
+            PrintWriter PW = new PrintWriter(new FileWriter(workingDir+"run"+proc+".bat"));
+            PW.println("python -u " + workingDir + "fetch"+proc+".py > " + workingDir+"out"+proc+".txt 2>&1");
+            PW.close();
+            processes[proc] = Runtime.getRuntime().exec(workingDir+"run"+proc+".bat");
+
+          }
+        }
+      }
+    }
+  }
+
+}

src/com/mrc/rasterise/era5_land/Extract.java

@@ -0,0 +1,45 @@
+package com.mrc.rasterise.era5_land;
+
+import java.io.File;
+import java.io.PrintWriter;
+import java.util.GregorianCalendar;
+
+public class Extract {
+
+  public static void extract_point(String temp_path, String precip_path, float[] lon, float[] lat, 
+                                   GregorianCalendar start, GregorianCalendar end,
+                                   String out_csv) throws Exception {
+
+    PrintWriter PW = new PrintWriter(new File(out_csv));
+    PW.println("lon,lat,year,month,day,t_min,t_mean,t_max,rh_min,rh_mean,rh_max,precip");
+    GregorianCalendar now = start;
+
+    int[] ls_to_era5_x = Coords.ls_to_era5_x();
+    int[] ls_to_era5_y = Coords.ls_to_era5_y();
+
+
+    while (!now.after(end)) {
+      int year = now.get(GregorianCalendar.YEAR);
+      int month = now.get(GregorianCalendar.MONTH) + 1;
+      int day = now.get(GregorianCalendar.DAY_OF_MONTH);
+      OneDay day_data = OneDay.get_day_data(temp_path, precip_path, year, month, day);
+      for (int p = 0; p < lon.length; p++) {
+        int era_x = ls_to_era5_x[Coords.lon_to_ls(lon[p])];
+        int era_y = ls_to_era5_y[Coords.lat_to_ls(lat[p])];
+        PW.println(lon[p] + "," + lat[p] + "," + year + "," + month + "," + day +"," +
+          day_data.tm_min[era_x][era_y] + "," +
+          day_data.tm_mean[era_x][era_y] + "," +
+          day_data.tm_max[era_x][era_y] + "," +
+          day_data.rh_min[era_x][era_y] + "," +
+          day_data.rh_mean[era_x][era_y] + "," +
+          day_data.rh_max[era_x][era_y] + "," +
+          day_data.precip[era_x][era_y]);
+      }
+      PW.flush();
+
+      now.add(GregorianCalendar.DATE,  1);
+    }
+    PW.close();
+  }
+
+}