WGS84 to Mercator conversion stability

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/FreeFormPointSet.java

@@ -251,8 +251,7 @@ public Envelope getWgsEnvelope () {
     public WebMercatorExtents getWebMercatorExtents () {
         final int DEFAULT_ZOOM = 9;
         Envelope wgsEnvelope = this.getWgsEnvelope();
-        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, DEFAULT_ZOOM);
-        return webMercatorExtents;
+        return WebMercatorExtents.forBufferedWgsEnvelope(wgsEnvelope, DEFAULT_ZOOM);
     }
 
     /** Construct a freeform point set containing one opportunity at each specified geographic coordinate. */

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

@@ -8,7 +8,6 @@
 import com.csvreader.CsvReader;
 import com.google.common.io.LittleEndianDataInputStream;
 import com.google.common.io.LittleEndianDataOutputStream;
-import org.apache.commons.math3.util.FastMath;
 import org.geotools.coverage.grid.GridCoverage2D;
 import org.geotools.coverage.grid.GridCoverageFactory;
 import org.geotools.coverage.grid.io.AbstractGridFormat;
@@ -64,14 +63,17 @@
 
 import static com.conveyal.gtfs.util.Util.human;
 import static com.conveyal.r5.common.GeometryUtils.checkWgsEnvelopeSize;
-import static com.google.common.base.Preconditions.checkArgument;
 import static com.google.common.base.Preconditions.checkState;
 import static java.lang.Double.parseDouble;
-import static org.apache.commons.math3.util.FastMath.atan;
-import static org.apache.commons.math3.util.FastMath.cos;
-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 java.lang.Math.PI;
+import static java.lang.Math.atan;
+import static java.lang.Math.cos;
+import static java.lang.Math.log;
+import static java.lang.Math.pow;
+import static java.lang.Math.sinh;
+import static java.lang.Math.tan;
+import static java.lang.Math.toDegrees;
+import static java.lang.Math.toRadians;
 
 /**
  * Class that represents a grid of opportunity counts in the spherical Mercator "projection" at a given zoom level.
@@ -122,7 +124,7 @@ public Grid (WebMercatorExtents extents) {
      * @param wgsEnvelope Envelope of grid, in absolute WGS84 lat/lon coordinates
      */
     public Grid (int zoom, Envelope wgsEnvelope) {
-        this(WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom));
+        this(WebMercatorExtents.forBufferedWgsEnvelope(wgsEnvelope, zoom));
     }
 
     public static class PixelWeight {
@@ -458,22 +460,30 @@ 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;
+        return xPixel / (pow(2, zoom) * 256) * 360 - 180;
     }
 
     /**
@@ -484,27 +494,36 @@ 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);
     }
 
     /**
-     * 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.
+     * Given an integer web Mercator pixel number, return the latitude in degrees of the north edge of that pixel at the
+     * given zoom level relative to the top (north) 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.
+     *
+     * TODO Profile this some time because we had versions using both FastMath and built-in Math functions.
+     * The difference should be less significant these days. Java now has a StrictMath and the normal Math is optimized.
      */
-    public static double pixelToCenterLat (int yPixel, int zoom) {
-        return pixelToLat(yPixel + 0.5, zoom);
+    public static double pixelToLat (double yPixel, int zoom) {
+        final double tile = yPixel / 256d;
+        return toDegrees(atan(sinh(PI - tile * PI * 2 / pow(2, zoom))));
     }
 
     /**
-     * Return the latitude of the north edge of any pixel at the given zoom level and y coordinate relative to the top
-     * edge of the world (assuming an integer pixel). Noninteger pixels will return locations within the pixel.
-     * We're using FastMath here, because the built-in math functions were taking a large amount of time in profiling.
+     * 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 pixelToLat (double yPixel, int zoom) {
-        return FastMath.toDegrees(atan(sinh(Math.PI - (yPixel / 256d) / Math.pow(2, zoom) * 2 * Math.PI)));
+    public static double pixelToCenterLat (int yPixel, int zoom) {
+        return pixelToLat(yPixel + 0.5, zoom);
     }
 
     /**
@@ -604,7 +623,7 @@ public static List<Grid> fromCsv(InputStreamProvider csvInputStreamProvider,
         reader.close();
 
         checkWgsEnvelopeSize(envelope, "CSV points");
-        WebMercatorExtents extents = WebMercatorExtents.forWgsEnvelope(envelope, zoom);
+        WebMercatorExtents extents = WebMercatorExtents.forBufferedWgsEnvelope(envelope, zoom);
         checkPixelCount(extents, numericColumns.size());
 
         if (progressListener != null) {
@@ -676,7 +695,7 @@ public static List<Grid> fromShapefile (File shapefile, int zoom, ProgressListen
         ShapefileReader reader = new ShapefileReader(shapefile);
         Envelope envelope = reader.wgs84Bounds();
         checkWgsEnvelopeSize(envelope, "Shapefile");
-        WebMercatorExtents extents = WebMercatorExtents.forWgsEnvelope(envelope, zoom);
+        WebMercatorExtents extents = WebMercatorExtents.forBufferedWgsEnvelope(envelope, zoom);
         List<String> numericAttributes = reader.numericAttributes();
         Set<String> uniqueNumericAttributes = new HashSet<>(numericAttributes);
         if (uniqueNumericAttributes.size() != numericAttributes.size()) {
@@ -753,14 +772,12 @@ 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
         Envelope wgsEnvelope = freeForm.getWgsEnvelope();
-        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom);
+        WebMercatorExtents webMercatorExtents = WebMercatorExtents.forBufferedWgsEnvelope(wgsEnvelope, zoom);
         Grid grid = new Grid(webMercatorExtents);
         grid.name = freeForm.name;
         for (int f = 0; f < freeForm.featureCount(); f++) {

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,94 @@ 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);
+    }
+
+    /**
+     * 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