Improve performance for OLCI, S2 MSI and MERIS

idepix-core/src/main/java/org/esa/snap/idepix/core/util/IdepixUtils.java

@@ -147,6 +147,15 @@ public static GeoPos getGeoPos(GeoCoding geoCoding, int x, int y) {
      * @return the azimuth difference [degree]
      */
     public static double computeAzimuthDifference(final double vaa, final double saa) {
+        final double delta = (720.0 + vaa - saa) % 360.0;
+        if (delta > 180.0) {
+            return 360.0 - delta;
+        } else {
+            return delta;
+        }
+    }
+
+    public static double computeAzimuthDifference1(final double vaa, final double saa) {
         return MathUtils.RTOD * Math.acos(Math.cos(MathUtils.DTOR * (vaa - saa)));
     }
 }

idepix-meris/src/main/java/org/esa/snap/idepix/meris/IdepixMerisSlopeAspectOrientationOp.java

@@ -154,9 +154,9 @@ public void computeTileStack(Map<Band, Tile> targetTiles, Rectangle targetRectan
                     final float vza = viewZenithAngleTile.getSampleFloat(x, y);
                     final float vaa = viewAzimuthAngleTile.getSampleFloat(x, y);
                     final float saa = sunAzimuthAngleTile.getSampleFloat(x, y);
-                    if (x == 2783 && y == 642) {
-                        System.out.println("x, y = " + x + ", " + y);  // small subset, shadow
-                    }
+//                    if (x == 2783 && y == 642) {
+//                        System.out.println("x, y = " + x + ", " + y);  // small subset, shadow
+//                    }
                     final float[] latitudeDataMacropixel =
                             SlopeAspectOrientationUtils.get3x3MacropixelData(latitudeTile, y, x);
                     final float[] longitudeDataMacropixel =

idepix-meris/src/main/java/org/esa/snap/idepix/meris/IdepixMerisUtils.java

@@ -97,24 +97,23 @@ static Product computeCloudTopPressureProduct(Product sourceProduct) {
     }
 
     static double computeApparentSaa(double sza, double saa, double oza, double oaa, double lat) {
-        final double szaRad = sza * MathUtils.DTOR;
-        final double ozaRad = oza * MathUtils.DTOR;
-
-        double deltaPhi;
-        if (oaa < 0.0) {
-            deltaPhi = 360.0 - Math.abs(oaa) - saa;
-        } else {
-            deltaPhi = saa - oaa;
-        }
-        final double deltaPhiRad = deltaPhi * MathUtils.DTOR;
-        final double numerator = Math.tan(szaRad) - Math.tan(ozaRad) * Math.cos(deltaPhiRad);
-        final double denominator = Math.sqrt(Math.tan(ozaRad) * Math.tan(ozaRad) + Math.tan(szaRad) * Math.tan(szaRad) -
-                2.0 * Math.tan(szaRad) * Math.tan(ozaRad) * Math.cos(deltaPhiRad));
-
+//        final double deltaPhi;
+//        if (oaa < 0.0) {
+//            deltaPhi = 360.0 - Math.abs(oaa) - saa;
+//        } else {
+//            deltaPhi = saa - oaa;
+//        }
+//        final double cosDeltaPhi = Math.cos(deltaPhi * MathUtils.DTOR);
+        final double tanSza = Math.tan(sza * MathUtils.DTOR);
+        final double tanOza = Math.tan(oza * MathUtils.DTOR);
+        final double cosDeltaPhi = Math.cos((saa - oaa) * MathUtils.DTOR);
+        final double numerator = tanSza - tanOza * cosDeltaPhi;
+        final double denominator = Math.sqrt(tanOza * tanOza + tanSza * tanSza -
+                2.0 * tanSza * tanOza * cosDeltaPhi);
         double delta = Math.acos(numerator / denominator);
-// Sun in the North (Southern hemisphere), change sign!
+        // Sun in the North (Southern hemisphere), change sign!
         if (saa > 270. || saa < 90){
-            delta = -1.0 * delta;
+            delta = -delta;
         }
         if (oaa < 0.0) {
             return saa - delta * MathUtils.RTOD;

idepix-meris/src/main/java/org/esa/snap/idepix/meris/IdepixMerisWaterClassificationOp.java

@@ -333,9 +333,27 @@ private double computeChiW(int x, int y, Tile winduTile, Tile windvTile, Tile sa
         return MathUtils.RTOD * (Math.acos(Math.cos(arg)));
     }
 
+    static double computeChiW1(double windU, double windV, double saa) {
+        final double phiw = azimuth(windU, windV);
+        /* and "scattering" angle */
+        final double arg = MathUtils.DTOR * (saa - phiw);
+        return MathUtils.RTOD * (Math.acos(Math.cos(arg)));
+    }
+
+    static double computeChiW2(double windU, double windV, double saa) {
+        final double phiw = MathUtils.RTOD * Math.atan2(windU, windV);
+        final double delta = (720.0 + saa - phiw) % 360.0;
+        if (delta > 180.0) {
+            return 360.0 - delta;
+        } else {
+            return delta;
+        }
+    }
+
     private double computeRhoGlint(int x, int y, Tile winduTile, Tile windvTile,
                                    Tile szaTile, Tile vzaTile, Tile saaTile, Tile vaaTile) {
-        final double chiw = computeChiW(x, y, winduTile, windvTile, saaTile);
+        //final double chiw = computeChiW(x, y, winduTile, windvTile, saaTile);
+        final double chiw = computeChiW2(winduTile.getSampleFloat(x, y), windvTile.getSampleFloat(x, y), saaTile.getSampleFloat(x, y));
         final float vaa = vaaTile.getSampleFloat(x, y);
         final float saa = saaTile.getSampleFloat(x, y);
         final double deltaAzimuth = (float) IdepixUtils.computeAzimuthDifference(vaa, saa);
@@ -390,7 +408,7 @@ private GeoPos getGeoPos(int x, int y) {
         return geoPos;
     }
 
-    private double azimuth(double x, double y) {
+    static double azimuth1(double x, double y) {
         if (y > 0.0) {
             // DPM #2.6.5.1.1-1
             return (MathUtils.RTOD * Math.atan(x / y));
@@ -402,6 +420,9 @@ private double azimuth(double x, double y) {
             return (x >= 0.0 ? 90.0 : 270.0);
         }
     }
+    static double azimuth(double y, double x) {
+        return MathUtils.RTOD * Math.atan2(y, x);
+    }
 
     public static class Spi extends OperatorSpi {
         public Spi() {

idepix-meris/src/test/java/org/esa/snap/idepix/meris/IdepixMerisWaterClassificationOpTest.java

@@ -0,0 +1,63 @@
+package org.esa.snap.idepix.meris;
+
+import junit.framework.TestCase;
+import org.esa.snap.core.util.math.MathUtils;
+import org.esa.snap.idepix.core.util.IdepixUtils;
+
+/**
+ * TODO add API doc
+ *
+ * @author Martin Boettcher
+ */
+public class IdepixMerisWaterClassificationOpTest extends TestCase {
+    public void testAtan() {
+        double ator = Math.PI / 180.0;
+        for (double y = -10.0; y <= 10.0; y += 1.0) {
+            for (double x = -10.0; y <= 10.0; y += 1.0) {
+                double v1 = IdepixMerisWaterClassificationOp.azimuth1(x * ator, y * ator);
+                double v2 = IdepixMerisWaterClassificationOp.azimuth(x * ator, y * ator);
+                assertEquals(x+","+y, v1, v2, 1.0e-5);
+            }
+        }
+    }
+
+    public void testChiw() {
+        double ator = Math.PI / 180.0;
+        for (double u = -180.0; u <= 180.0; u += 45.0) {
+            for (double v = -180.0; v <= 180.0; v += 45.0) {
+                for (double s = -180.0; s <= 180.0; s += 15.0) {
+                    double w1 = IdepixMerisWaterClassificationOp.computeChiW1(u * ator, v * ator, s);
+                    double w2 = IdepixMerisWaterClassificationOp.computeChiW2(u * ator, v * ator, s);
+                    assertEquals(u + "," + v + "," + s, w1, w2, 1.0e-5);
+                }
+            }
+        }
+    }
+
+    public void testCos() {
+        for (double oaa = -180.0; oaa <= 360.0; oaa += 45.0) {
+            for (double saa = -180.0; saa <= 360.0; saa += 45.0) {
+                final double deltaPhi;
+                if (oaa < 0.0) {
+                    deltaPhi = 360.0 - Math.abs(oaa) - saa;
+                } else {
+                    deltaPhi = saa - oaa;
+                }
+                double v1 = Math.cos(MathUtils.DTOR * deltaPhi);
+                double v2 = Math.cos(MathUtils.DTOR * (saa - oaa));
+                assertEquals(oaa + "," + saa, v1, v2, 1e-5);
+            }
+        }
+    }
+
+    public void testAzimuthDifference() {
+        for (double oaa = -180.0; oaa <= 360.0; oaa += 45.0) {
+            for (double saa = -180.0; saa <= 360.0; saa += 45.0) {
+                double v1 = IdepixUtils.computeAzimuthDifference1(oaa, saa);
+                double v2 = IdepixUtils.computeAzimuthDifference(oaa, saa);
+                assertEquals(oaa + "," + saa, v1, v2, 1e-5);
+            }
+        }
+    }
+
+}

idepix-olci/pom.xml

@@ -26,7 +26,7 @@
     </parent>
 
     <artifactId>idepix-olci</artifactId>
-    <version>9.0.2-SNAPSHOT</version>
+    <version>9.0.3-SNAPSHOT</version>
 
     <packaging>nbm</packaging>
 

idepix-olci/src/main/java/org/esa/snap/idepix/olci/CtpOp.java

@@ -139,8 +139,11 @@ public void doExecute(ProgressMonitor pm) throws OperatorException {
 
     @Override
     public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) throws OperatorException {
+        if (!"ctp".equals(targetBand.getName())) {
+            throw new OperatorException("Unexpected target band name: '" + targetBand.getName() + "' - exiting.");
+        }
+
         final Rectangle targetRectangle = targetTile.getRectangle();
-        final String targetBandName = targetBand.getName();
 
         final Tile szaTile = getSourceTile(szaBand, targetRectangle);
         final Tile ozaTile = getSourceTile(ozaBand, targetRectangle);
@@ -154,6 +157,8 @@ public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) th
 
         final Tile l1FlagsTile = getSourceTile(sourceProduct.getRasterDataNode("quality_flags"), targetRectangle);
 
+        float[][] nnInputs = new float[targetRectangle.height * targetRectangle.width][];
+        float[] dummyNnInput = new float[7];
         for (int y = targetRectangle.y; y < targetRectangle.y + targetRectangle.height; y++) {
             checkForCancellation();
             for (int x = targetRectangle.x; x < targetRectangle.x + targetRectangle.width; x++) {
@@ -180,22 +185,33 @@ public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) th
                     final float tra15 = tra15Tile.getSampleFloat(x, y);
                     final float mLogTra15 = (float) -Math.log(tra15);
 
-                    float[] nnInput = new float[]{cosSza, cosOza, aziDiff, refl12, mLogTra13, mLogTra14, mLogTra15};
-                    final float[][] nnResult = nnCalculator.calculate(nnInput);
-                    final float ctp = TensorflowNNCalculator.convertNNResultToCtp(nnResult[0][0]);
+                    float[] nnInput = new float[] {cosSza, cosOza, aziDiff, refl12, mLogTra13, mLogTra14, mLogTra15};
+                    //final float[][] nnResult = nnCalculator.calculate(nnInput);
+                    //final float ctp = TensorflowNNCalculator.convertNNResultToCtp(nnResult[0][0]);
+                    //targetTile.setSample(x, y, ctp);
+                    nnInputs[(y-targetRectangle.y) * targetRectangle.width + (x-targetRectangle.x)] = nnInput;
+                } else {
+                    //targetTile.setSample(x, y, Float.NaN);
+                    nnInputs[(y-targetRectangle.y) * targetRectangle.width + (x-targetRectangle.x)] = dummyNnInput;
+                }
+            }
+        }
+
+        // call tensorflow once with the complete tile stack
+        final float[][] nnResult = nnCalculator.calculate1(nnInputs);
 
-                    if (targetBandName.equals("ctp")) {
-                        targetTile.setSample(x, y, ctp);
-                    } else {
-                        throw new OperatorException("Unexpected target band name: '" +
-                                                            targetBandName + "' - exiting.");
-                    }
+        // convert output of tf into ctp and set value into target tile
+        for (int y = targetRectangle.y; y < targetRectangle.y + targetRectangle.height; y++) {
+            checkForCancellation();
+            for (int x = targetRectangle.x; x < targetRectangle.x + targetRectangle.width; x++) {
+                final boolean pixelIsValid = !l1FlagsTile.getSampleBit(x, y, IdepixOlciConstants.L1_F_INVALID);
+                if (pixelIsValid) {
+                    targetTile.setSample(x, y, TensorflowNNCalculator.convertNNResultToCtp(nnResult[(y-targetRectangle.y) * targetRectangle.width + (x-targetRectangle.x)][0]));
                 } else {
                     targetTile.setSample(x, y, Float.NaN);
                 }
             }
         }
-
     }
 
     private void preProcess() {

idepix-olci/src/main/java/org/esa/snap/idepix/olci/TensorflowNNCalculator.java

@@ -40,9 +40,11 @@ class TensorflowNNCalculator {
         try {
             TensorFlow.version(); // triggers init of TensorFlow
         } catch (LinkageError e) {
-            throw new IllegalStateException("TensorFlow could not be initialised. " +
-                                                    "Make sure that your CPU supports 64Bit and AVX instruction set " +
-                                                    "(Are you using a VM?) and that you have installed the Microsoft Visual C++ 2015 Redistributable when you are on windows.", e);
+            throw new IllegalStateException(
+                    "TensorFlow could not be initialised. " +
+                            "Make sure that your CPU supports 64Bit and AVX instruction set " +
+                            "(Are you using a VM?) and that you have installed the Microsoft Visual C++ 2015 " +
+                            "Redistributable when you are on windows.", e);
         }
 
         this.transformMethod = transformMethod;
@@ -181,7 +183,6 @@ private void loadModel() throws Exception {
      * Requires that loadModel() is run once before.
      *
      * @param nnInput - input vector for neural net
-     *
      * @return float[][] - the converted result array
      */
     float[][] calculate(float[] nnInput) {
@@ -210,5 +211,42 @@ float[][] calculate(float[] nnInput) {
         }
     }
 
+    /**
+     * Applies NN to vector of pixel band stacks and returns converted array.
+     * Functional implementation of setNnTensorInput(.) plus getNNResult().
+     * Makes sure the Tensors are closed after use.
+     * Requires that loadModel() is run once before.
+     *
+     * @param nnInput - image vector of band vectors, band vectors are input for neural net
+     * @return float[][] - image vector of output band vector (length 1)
+     */
+    float[][] calculate1(float[][] nnInput) {
+        if (transformMethod.equals("sqrt")) {
+            for (int i = 0; i < nnInput.length; i++) {
+                for (int j = 0; j < nnInput[i].length; j++) {
+                    nnInput[i][j] = (float) Math.sqrt(nnInput[i][j]);
+                }
+            }
+        } else if (transformMethod.equals("log")) {
+            for (int i = 0; i < nnInput.length; i++) {
+                for (int j = 0; j < nnInput[i].length; j++) {
+                    nnInput[i][j] = (float) Math.log10(nnInput[i][j]);
+                }
+            }
+        }
+        final Session.Runner runner = model.session().runner();
+        try (
+                Tensor<?> inputTensor = Tensor.create(nnInput);
+                Tensor<?> outputTensor = runner.feed(firstNodeName, inputTensor).fetch(lastNodeName).run().get(0)
+        ) {
+            long[] ts = outputTensor.shape();
+            int numPixels = (int) ts[0];
+            int numOutputVars = (int) ts[1];
+            float[][] m = new float[numPixels][numOutputVars];
+            outputTensor.copyTo(m);
+            return m;
+        }
+    }
+
 
 }

idepix-s2msi/pom.xml

@@ -26,7 +26,7 @@
     </parent>
 
     <artifactId>idepix-s2msi</artifactId>
-    <version>9.0.2-SNAPSHOT</version>
+    <version>9.0.3-SNAPSHOT</version>
 
     <packaging>nbm</packaging>
 

idepix-s2msi/src/main/java/org/esa/snap/idepix/s2msi/S2IdepixOp.java

@@ -11,6 +11,7 @@
 import org.esa.snap.core.gpf.annotations.SourceProduct;
 import org.esa.snap.core.gpf.annotations.TargetProduct;
 import org.esa.snap.dem.gpf.AddElevationOp;
+import org.esa.snap.idepix.s2msi.operators.S2IdepixCloudPostProcess2Op;
 import org.esa.snap.idepix.s2msi.operators.S2IdepixCloudPostProcessOp;
 import org.esa.snap.idepix.s2msi.util.AlgorithmSelector;
 import org.esa.snap.idepix.s2msi.util.S2IdepixConstants;
@@ -108,6 +109,10 @@ private void processSentinel2() {
         int cacheSize = Integer.parseInt(System.getProperty(S2IdepixUtils.TILECACHE_PROPERTY, "1600"));
         s2ClassifProduct = S2IdepixUtils.computeTileCacheProduct(s2ClassifProduct, cacheSize);
 
+        // breakpoint output to generate input for cloud post-processing
+        //targetProduct = s2ClassifProduct;
+        //if (true) return;
+
         // Post Cloud Classification: cloud shadow, cloud buffer, mountain shadow
         Product postProcessingProduct = computePostProcessProduct(sourceProduct, s2ClassifProduct);
 
@@ -146,7 +151,10 @@ private Product computePostProcessProduct(Product l1cProduct, Product classifica
         Map<String, Object> paramsBuffer = new HashMap<>();
         paramsBuffer.put("cloudBufferWidth", cloudBufferWidth);
         paramsBuffer.put("computeCloudBufferForCloudAmbiguous", computeCloudBufferForCloudAmbiguous);
-        Product cloudBufferProduct = GPF.createProduct(OperatorSpi.getOperatorAlias(S2IdepixCloudPostProcessOp.class),
+        //Product cloudBufferProduct = GPF.createProduct(OperatorSpi.getOperatorAlias(S2IdepixCloudPostProcessOp.class),
+        //                                               paramsBuffer, input);
+        paramsBuffer.put("computeCloudBuffer", computeCloudBuffer);
+        Product cloudBufferProduct = GPF.createProduct(OperatorSpi.getOperatorAlias(S2IdepixCloudPostProcess2Op.class),
                                                        paramsBuffer, input);
 
         int cacheSize = Integer.parseInt(System.getProperty(S2IdepixUtils.TILECACHE_PROPERTY, "1600")) / 5;

idepix-s2msi/src/main/java/org/esa/snap/idepix/s2msi/UrbanCloudDistinction.java

@@ -163,17 +163,21 @@ public void correctCloudFlag(int x, int y, Tile classifFlagTile, Tile targetFlag
 
     }
 
+    public float getCdiValue(int x, int y) {
+        return cdiBand.getSampleFloat(x, y);
+    }
+
 
     /**
      * Adds debug band to the target product.
      */
-    public void addDebugBandsToTargetProduct() {
-        ucd_oldCloud_debug = s2ClassifProduct.addBand("__ucd_oldCloud__", ProductData.TYPE_INT8);
+    public void addDebugBandsToTargetProduct(Product targetProduct) {
+        ucd_oldCloud_debug = targetProduct.addBand("__ucd_oldCloud__", ProductData.TYPE_INT8);
         ucd_oldCloud_debug.ensureRasterData();
-        ucd_cdi_debug = s2ClassifProduct.addBand("__ucd_cdi__", ProductData.TYPE_FLOAT32);
+        ucd_cdi_debug = targetProduct.addBand("__ucd_cdi__", ProductData.TYPE_FLOAT32);
         ucd_cdi_debug.setNoDataValue(Float.NaN);
         ucd_cdi_debug.ensureRasterData();
-        ucd_cloudMean11_debug = s2ClassifProduct.addBand("__ucd_cloudMean11__", ProductData.TYPE_FLOAT32);
+        ucd_cloudMean11_debug = targetProduct.addBand("__ucd_cloudMean11__", ProductData.TYPE_FLOAT32);
         ucd_cloudMean11_debug.setNoDataValue(Float.NaN);
         ucd_cloudMean11_debug.ensureRasterData();
         debugBandsEnabled = true;