diff --git a/src/main/java/com/conveyal/analysis/models/AnalysisRequest.java b/src/main/java/com/conveyal/analysis/models/AnalysisRequest.java index b8859d9d5..80d5e507a 100644 --- a/src/main/java/com/conveyal/analysis/models/AnalysisRequest.java +++ b/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; diff --git a/src/main/java/com/conveyal/r5/analyst/Grid.java b/src/main/java/com/conveyal/r5/analyst/Grid.java index 60bb526c6..a47021b0f 100644 --- a/src/main/java/com/conveyal/r5/analyst/Grid.java +++ b/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 diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java index 0d75a79ca..cff116c00 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java +++ b/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 diff --git a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java new file mode 100644 index 000000000..e0600cb6e --- /dev/null +++ b/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); + } + +} \ No newline at end of file