Revise WebMercatorExtents factory methods

Handle edge cases where the WGS84 envelope lies on Mercator pixel edges.
Handle numerical instability in repeated coordinate system conversions.
Add tests for repeated coordinate conversion.

src/main/java/com/conveyal/r5/analyst/Grid.java

@@ -458,19 +458,25 @@ public int featureCount() {
         return extents.width * extents.height;
     }
 
-    /* functions below from http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Mathematics */
+    /* Functions below derived from http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Mathematics */
+
+    /** Like lonToPixel but returns fractional values for positions within a pixel instead of integer pixel numbers. */
+    public static double lonToFractionalPixel (double lon, int zoom) {
+        return (lon + 180) / 360 * Math.pow(2, zoom) * 256;
+    }
 
     /**
-     * Return the absolute (world) x pixel number of all pixels the given line of longitude falls within at the given
-     * zoom level.
+     * Return the absolute (world) x pixel number of all pixels containing the given line of longitude
+     * at the given zoom level.
      */
     public static int lonToPixel (double lon, int zoom) {
-        return (int) ((lon + 180) / 360 * Math.pow(2, zoom) * 256);
+        return (int) lonToFractionalPixel(lon, zoom);
     }
 
     /**
      * Return the longitude of the west edge of any pixel at the given zoom level and x pixel number measured from the
-     * west edge of the world (assuming an integer pixel). Noninteger pixels will return locations within that pixel.
+     * west edge of the world (assuming an integer x pixel number). Noninteger x pixel coordinates will return
+     * WGS84 locations within that pixel.
      */
     public static double pixelToLon (double xPixel, int zoom) {
         return xPixel / (Math.pow(2, zoom) * 256) * 360 - 180;
@@ -484,18 +490,15 @@ public static double pixelToCenterLon (int xPixel, int zoom) {
         return pixelToLon(xPixel + 0.5, zoom);
     }
 
-    /** Return the absolute (world) y pixel number of all pixels the given line of latitude falls within. */
-    public static int latToPixel (double lat, int zoom) {
+    /** Like latToPixel but returns fractional values for positions within a pixel instead of integer pixel numbers. */
+    public static double latToFractionalPixel (double lat, int zoom) {
         double latRad = FastMath.toRadians(lat);
-        return (int) ((1 - log(tan(latRad) + 1 / cos(latRad)) / Math.PI) * Math.pow(2, zoom - 1) * 256);
+        return (1 - log(tan(latRad) + 1 / cos(latRad)) / Math.PI) * Math.pow(2, zoom - 1) * 256;
     }
 
-    /**
-     * Return the latitude of the center of all pixels at the given zoom level and absolute (world) y pixel number
-     * measured southward from the north edge of the world.
-     */
-    public static double pixelToCenterLat (int yPixel, int zoom) {
-        return pixelToLat(yPixel + 0.5, zoom);
+    /** Return the absolute (world) y pixel number of all pixels the given line of latitude falls within. */
+    public static int latToPixel (double lat, int zoom) {
+        return (int) latToFractionalPixel(lat, zoom);
     }
 
     /**
@@ -507,6 +510,14 @@ public static double pixelToLat (double yPixel, int zoom) {
         return FastMath.toDegrees(atan(sinh(Math.PI - (yPixel / 256d) / Math.pow(2, zoom) * 2 * Math.PI)));
     }
 
+    /**
+     * Return the latitude of the center of all pixels at the given zoom level and absolute (world) y pixel number
+     * measured southward from the north edge of the world.
+     */
+    public static double pixelToCenterLat (int yPixel, int zoom) {
+        return pixelToLat(yPixel + 0.5, zoom);
+    }
+
     /**
      * Given a pixel's local grid coordinates within the supplied WebMercatorExtents, return a closed
      * polygon of that pixel's outline in WGS84 global geographic coordinates.

src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java

@@ -10,8 +10,7 @@
 
 import java.util.Arrays;
 
-import static com.conveyal.r5.analyst.Grid.latToPixel;
-import static com.conveyal.r5.analyst.Grid.lonToPixel;
+import static com.conveyal.r5.analyst.Grid.*;
 import static com.google.common.base.Preconditions.checkArgument;
 import static com.google.common.base.Preconditions.checkElementIndex;
 import static com.google.common.base.Preconditions.checkNotNull;
@@ -114,23 +113,91 @@ public WebMercatorExtents expandToInclude (WebMercatorExtents other) {
         return new WebMercatorExtents(outWest, outNorth, outWidth, outHeight, this.zoom);
     }
 
+    /**
+     * Construct a WebMercatorExtents that contains all the points inside the supplied WGS84 envelope.
+     * The logic here is designed to handle the case where the edges of the WGS84 envelope exactly coincide with the
+     * edges of the web Mercator grid. For example, if the right edge is exactly on the border of pixels with x
+     * number 3 and 4, and the left edge is exactly on the border of pixels with x number 2 and 3, the resulting
+     * WebMercatorExtents should only include pixels with x number 3. This case arises when the WGS84 envelope was
+     * itself derived from a block of web Mercator pixels.
+     *
+     * Note that this does not handle numerical instability or imprecision in repeated floating point calculations on
+     * envelopes, where WGS84 coordinates are intended to yield exact integer web Mercator coordinates but do not.
+     * Such imprecision is handled by wrapper factory methods (by shrinking or growing the envelope by some epsilon).
+     * The real solution would be to always specify origin grids in integer form, but our current API does not allow this.
+     */
     public static WebMercatorExtents forWgsEnvelope (Envelope wgsEnvelope, int zoom) {
-        /*
-          The grid extent is computed from the points. If the cell number for the right edge of the grid is rounded
-          down, some points could fall outside the grid. `latToPixel` and `lonToPixel` naturally truncate down, which is
-          the correct behavior for binning points into cells but means the grid is (almost) always 1 row too
-          narrow/short, so we add 1 to the height and width when a grid is created in this manner. The exception is
-          when the envelope edge lies exactly on a pixel boundary. For this reason we should probably not produce WGS
-          Envelopes that exactly align with pixel edges, but they should instead surround the points at pixel centers.
-          Note also that web Mercator coordinates increase from north to south, so minimum latitude is maximum y.
-          TODO maybe use this method when constructing Grids. Grid (int zoom, Envelope envelope)
-         */
+        // Conversion of WGS84 to web Mercator pixels truncates intra-pixel coordinates toward the origin (northwest).
+        // Note that web Mercator coordinates increase from north to south, so maximum latitude is minimum Mercator y.
         int north = latToPixel(wgsEnvelope.getMaxY(), zoom);
         int west = lonToPixel(wgsEnvelope.getMinX(), zoom);
-        int height = (latToPixel(wgsEnvelope.getMinY(), zoom) - north) + 1; // minimum height is 1
-        int width = (lonToPixel(wgsEnvelope.getMaxX(), zoom) - west) + 1; // minimum width is 1
-        WebMercatorExtents webMercatorExtents = new WebMercatorExtents(west, north, width, height, zoom);
-        return webMercatorExtents;
+        // Find width and height in whole pixels, handling the case where the right envelope edge is on a pixel edge.
+        int height = (int) Math.ceil(latToFractionalPixel(wgsEnvelope.getMinY(), zoom) - north);
+        int width = (int) Math.ceil(lonToFractionalPixel(wgsEnvelope.getMaxX(), zoom) - west);
+        // These extents are constructed to contain some objects, and they have integer sizes (whole Mercator pixels).
+        // Therefore they should always have positive nonzero integer width and height.
+        checkState(width >= 1, "Web Mercator extents should always have a width of one or more.");
+        checkState(height >= 1, "Web Mercator extents should always have a height of one or more.");
+        return new WebMercatorExtents(west, north, width, height, zoom);
+    }
+
+    /**
+     * For function definitions below to work, this epsilon must be (significantly) less than half the height or
+     * width of a web Mercator pixel at the highest zoom level and highest latitude allowed by the system.
+     * Current value is about 1 meter north-south, but east-west distance is higher at higher latitudes.
+     */
+    private static final double WGS_EPSILON = 0.00001;
+
+    /**
+     * Wrapper to forWgsEnvelope where the envelope is shrunk uniformly by only a meter or two in each direction. This
+     * helps handle lack of numerical precision where floating point math is intended to yield integers and envelopes
+     * that exactly line up with Mercator grid edges (though the latter is also handled by the main constructor).
+     * A left or top edge coordinate that ends up the tiniest bit below the intended integer will be truncated all the
+     * way to the next lower pixel number. A bottom or right edge coordinate that ends up a tiny bit above the intended
+     * integer will tack a whole row of pixels onto the grid. Neither of these are horrible in isolation (grid is one
+     * pixel too big) but when applied repeatedly, as when basing one analysis on another in a chain, the grid will grow
+     * larger and larger, yielding incomparable result grids.
+     *
+     * I originally attempted to do this with lat/lonToPixel functions that would snap to pixel edges, but the snapping
+     * needs to pull in different directions on different edges of the envelope. Just transforming the envelope is simpler.
+     *
+     * Trimming is not appropriate in every use case. If you have an envelope that's a tight fit around a set of points
+     * (opportunities) and the left edge of that envelope is within the trimming distance of a web Mercator pixel edge,
+     * those opportunities will be outside the resulting WebMercatorExtents. When preparing grids to contain sets of
+     * points, it is probably better to deal with numerical imprecision by expanding the envelope slightly instead of
+     * contracting it.
+     */
+    public static WebMercatorExtents forTrimmedWgsEnvelope (Envelope wgsEnvelope, int zoom) {
+        checkArgument(wgsEnvelope.getWidth() > 3 * WGS_EPSILON, "Envelope is too narrow.");
+        checkArgument(wgsEnvelope.getHeight() > 3 * WGS_EPSILON, "Envelope is too short.");
+        // Shrink a protective copy of the envelope slightly. Note the negative sign, this is shrinking the envelope.
+        wgsEnvelope = wgsEnvelope.copy();
+        wgsEnvelope.expandBy(-WGS_EPSILON);
+        return WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom);
+    }
+
+    /**
+     * The opposite of forTrimmedWgsEnvelope: makes the envelope a tiny bit bigger before constructing the extents.
+     * This helps deal with numerical imprecision where we want to be sure all points within a supplied envelope will
+     * fall inside cells of the resulting web Mercator grid.
+     */
+    public static WebMercatorExtents forBufferedWgsEnvelope (Envelope wgsEnvelope, int zoom) {
+        // Expand a protective copy of the envelope slightly.
+        wgsEnvelope = wgsEnvelope.copy();
+        wgsEnvelope.expandBy(WGS_EPSILON);
+        return WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom);
+    }
+
+    /**
+     * Produces a new Envelope in WGS84 coordinates that tightly encloses the pixels of this WebMercatorExtents.
+     * The edges of that Envelope will run exactly down the borders between neighboring web Mercator pixels.
+     */
+    public Envelope toWgsEnvelope () {
+        double westLon = pixelToLon(west, zoom);
+        double northLat = pixelToLat(north, zoom);
+        double eastLon = pixelToLon(west + width, zoom);
+        double southLat = pixelToLat(north + height, zoom);
+        return new Envelope(westLon, eastLon, southLat, northLat);
     }
 
     @Override

src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java

@@ -150,6 +150,10 @@ public TIntList getPointsInEnvelope (Envelope envelopeFixedDegrees) {
     }
 
     // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Mathematics
+    // FIXME Grid.java has seemingly identical definitions of these same mercator-to-wgs methods.
+    //   These methods should just be wrappers that pass in this WebMercatorGridPointSet's zoom level.
+    //   Grid also has methods that distinguish between pixel corner and center. Those corner methods should be renamed.
+    //   WebMercatorGridPointSet, Grid, and WebMercatorExtents are in the same package and should share methods.
 
     /**
      * Return the x pixel number containing the given longitude at this grid's zoom level, relative to the left edge

src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java

@@ -0,0 +1,65 @@
+package com.conveyal.r5.analyst;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+import org.locationtech.jts.geom.Envelope;
+
+import static org.junit.jupiter.api.Assertions.*;
+
+class WebMercatorExtentsTest {
+
+    /**
+     * Normal usage of Conveyal Analysis involves a feedback loop between Web Mercator extents and WGS84 extents.
+     * This is because the API endpoint to create a new regional analysis requires extents in WGS84, but these are
+     * always converted into Web Mercator to run the actual analysis. If someone wants to run another analysis with
+     * the same dimensions, the UI must extract them from the database record or results of the previous analysis and
+     * pass back through the WGS84 coordinate system.
+     *
+     * We originally made an ill-advised assumption that most WGS84 extents would not fall exactly on Mercator pixel
+     * boundaries. This is sort of true for locations drawn from nature: the vast majority of real numbers are not
+     * integers. But once we discretize these into Mercator pixels/cells and feed those values back into the conversion,
+     * of course they are integers! In addition, there may be numerical instability in the conversion such that we
+     * sometimes get integers and other times coordinates falling ever so slightly on either side of the pixel boundary.
+     *
+     * This test verifies that such repeated conversions are stable and do not yield an ever-increasing envelope size.
+     * Ideally we should also test that this is the case when the non-integer WGS84 values are truncated by a few
+     * digits (as when they are stored in a database or serialized over the wire).
+     */
+    @Test
+    void wgsMercatorStability () {
+        Envelope wgsEnvelope = new Envelope();
+        wgsEnvelope.expandToInclude(10.22222222, 45.111111);
+        wgsEnvelope.expandBy(0.1);
+        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, 9);
+
+        for (int i = 0; i < 10; i++) {
+            Assertions.assertTrue(webMercatorExtents.width > 20);
+            Assertions.assertTrue(webMercatorExtents.height > 20);
+            wgsEnvelope = webMercatorExtents.toWgsEnvelope();
+            // Note the use of the trimmed envelope factory function, which should be used in the API.
+            WebMercatorExtents webMercatorExtents2 = WebMercatorExtents.forTrimmedWgsEnvelope(wgsEnvelope, 9);
+            Assertions.assertEquals(webMercatorExtents2, webMercatorExtents);
+            webMercatorExtents = webMercatorExtents2;
+        }
+    }
+
+    /**
+     * Check that a zero-size envelope (around a single point for example) will yield an extents object containing
+     * one cell (rather than zero cells). Also check an envelope with a tiny nonzero envelope away from cell edges.
+     */
+    @Test
+    void singleCellExtents () {
+        Envelope wgsEnvelope = new Envelope();
+
+        wgsEnvelope.expandToInclude(10.22222222, 10.32222222);
+        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, 9);
+        Assertions.assertEquals(1, webMercatorExtents.width);
+        Assertions.assertEquals(1, webMercatorExtents.height);
+
+        wgsEnvelope.expandBy(0.00001);
+        webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, 9);
+        Assertions.assertEquals(1, webMercatorExtents.width);
+        Assertions.assertEquals(1, webMercatorExtents.height);
+    }
+
+}