Make repeated WGS84-Mercator conversion stable

src/main/java/com/conveyal/analysis/models/AnalysisRequest.java

@@ -216,9 +216,7 @@ public void populateTask (AnalysisWorkerTask task, UserPermissions userPermissio
         task.maxFare = maxFare;
         task.inRoutingFareCalculator = inRoutingFareCalculator;
 
-        // TODO define class with static factory function WebMercatorGridBounds.fromLatLonBounds().
-        //      Also include getIndex(x, y), getX(index), getY(index), totalTasks()
-        WebMercatorExtents extents = WebMercatorExtents.forWgsEnvelope(bounds.envelope(), zoom);
+        WebMercatorExtents extents = WebMercatorExtents.forTrimmedWgsEnvelope(bounds.envelope(), zoom);
         task.height = extents.height;
         task.north = extents.north;
         task.west = extents.west;

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

@@ -72,6 +72,7 @@
 import static org.apache.commons.math3.util.FastMath.log;
 import static org.apache.commons.math3.util.FastMath.sinh;
 import static org.apache.commons.math3.util.FastMath.tan;
+import static org.apache.commons.math3.util.FastMath.toRadians;
 
 /**
  * Class that represents a grid of opportunity counts in the spherical Mercator "projection" at a given zoom level.
@@ -458,19 +459,27 @@ 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) {
+        // Factor of 256 yields a pixel value rather than the tile number.
+        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.
+     * Given an integer web Mercator pixel number, return the longitude in degrees of the west edge of that pixel at
+     * the given zoom level measured from the west edge of the world (not relative to this grid or to any particular
+     * tile). Given a non-integer web Mercator pixel number, return WGS84 locations within that pixel.
+     * Naming should somehow be revised to clarify that it doesn't return the center of the pixel.
      */
     public static double pixelToLon (double xPixel, int zoom) {
         return xPixel / (Math.pow(2, zoom) * 256) * 360 - 180;
@@ -484,10 +493,15 @@ public static double pixelToCenterLon (int xPixel, int zoom) {
         return pixelToLon(xPixel + 0.5, 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) {
+        final double latRad = toRadians(lat);
+        return (1 - log(tan(latRad) + 1 / cos(latRad)) / Math.PI) * Math.pow(2, zoom - 1) * 256;
+    }
+
     /** 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) {
-        double latRad = FastMath.toRadians(lat);
-        return (int) ((1 - log(tan(latRad) + 1 / cos(latRad)) / Math.PI) * Math.pow(2, zoom - 1) * 256);
+        return (int) latToFractionalPixel(lat, zoom);
     }
 
     /**
@@ -753,9 +767,7 @@ public double getOpportunityCount (int i) {
     }
 
     /**
-     * Rasterize a FreeFormFointSet into a Grid.
-     * Currently intended only for UI display of FreeFormPointSets, or possibly for previews of accessibility results
-     * during
+     * Rasterize a FreeFormFointSet into a Grid. Currently intended only for UI display of FreeFormPointSets.
      */
     public static Grid fromFreeForm (FreeFormPointSet freeForm, int zoom) {
         // TODO make and us a strongly typed WgsEnvelope class here and in various places

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

@@ -10,8 +10,12 @@
 
 import java.util.Arrays;
 
+import static com.conveyal.r5.analyst.Grid.latToFractionalPixel;
 import static com.conveyal.r5.analyst.Grid.latToPixel;
+import static com.conveyal.r5.analyst.Grid.lonToFractionalPixel;
 import static com.conveyal.r5.analyst.Grid.lonToPixel;
+import static com.conveyal.r5.analyst.Grid.pixelToLat;
+import static com.conveyal.r5.analyst.Grid.pixelToLon;
 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 +118,82 @@ 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 trimmed 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).
+     *
+     * Without this trimming, when receiving an envelope that's supposed to lie on integer web Mercator pixel
+     * boundaries, a left or top edge coordinate may end up the tiniest bit below the intended integer and be
+     * truncated all the way down 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 (the grid is just one pixel bigger than absolutely necessary) but when applied repeatedly, as when
+     * basing one analysis on another in a chain, the grid will grow larger and larger, yielding mismatched result grids.
+     *
+     * I originally attempted to handle 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
+     * by trimming it slightly is simpler and has the intended effect.
+     *
+     * 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 intended to contain
+     * a set of points, it is probably better to deal with numerical imprecision by expanding the envelope slightly
+     * instead of contracting it. A separate method is supplied for this purpose.
+     */
+    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);
+    }
+
+    /**
+     * 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/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java

@@ -0,0 +1,71 @@
+package com.conveyal.r5.analyst;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.ValueSource;
+import org.locationtech.jts.geom.Envelope;
+
+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).
+     */
+    @ParameterizedTest
+    @ValueSource(ints = {9, 10, 11, 12})
+    void wgsMercatorStability (int zoom) {
+        Envelope wgsEnvelope = new Envelope();
+        wgsEnvelope.expandToInclude(10.22222222, 45.111111);
+        wgsEnvelope.expandBy(0.1);
+        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom);
+
+        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, zoom);
+            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.
+     */
+    @ParameterizedTest
+    @ValueSource(ints = {9, 10, 11, 12})
+    void singleCellExtents (int zoom) {
+        Envelope wgsEnvelope = new Envelope();
+
+        wgsEnvelope.expandToInclude(10.22222222, 10.32222222);
+        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom);
+        Assertions.assertEquals(1, webMercatorExtents.width);
+        Assertions.assertEquals(1, webMercatorExtents.height);
+
+        wgsEnvelope.expandBy(0.00001);
+        webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom);
+        Assertions.assertEquals(1, webMercatorExtents.width);
+        Assertions.assertEquals(1, webMercatorExtents.height);
+
+        // Try taking these single-pixel extents through WGS84 for good measure.
+        WebMercatorExtents wme2 = WebMercatorExtents.forTrimmedWgsEnvelope(webMercatorExtents.toWgsEnvelope(), zoom);
+        Assertions.assertEquals(webMercatorExtents, wme2);
+    }
+
+}