From 519bf55193e1d623bcec57dec6cc546e41733abc Mon Sep 17 00:00:00 2001 From: Andrew Byrd Date: Mon, 16 Oct 2023 23:26:56 +0200 Subject: [PATCH 1/6] 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. --- .../java/com/conveyal/r5/analyst/Grid.java | 39 +++++--- .../r5/analyst/WebMercatorExtents.java | 99 ++++++++++++++++--- .../r5/analyst/WebMercatorGridPointSet.java | 4 + .../r5/analyst/WebMercatorExtentsTest.java | 65 ++++++++++++ 4 files changed, 177 insertions(+), 30 deletions(-) create mode 100644 src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java diff --git a/src/main/java/com/conveyal/r5/analyst/Grid.java b/src/main/java/com/conveyal/r5/analyst/Grid.java index 60bb526c6..628ec31ec 100644 --- a/src/main/java/com/conveyal/r5/analyst/Grid.java +++ b/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. diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java index 0d75a79ca..b17a181c0 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java +++ b/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 diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java index b6b614df7..46607a470 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java +++ b/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 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..5d0744556 --- /dev/null +++ b/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); + } + +} \ No newline at end of file From 5ae5fce5357e1797726da3387bd47ebe1b5998bb Mon Sep 17 00:00:00 2001 From: Andrew Byrd Date: Mon, 16 Oct 2023 23:37:40 +0200 Subject: [PATCH 2/6] Use new trimmed/buffered factory methods Trimming when handling possibly recycled WGS84 bounds. Buffering when supplying tight fit opportunity bounding boxes. --- .../conveyal/analysis/models/AnalysisRequest.java | 4 +--- .../com/conveyal/r5/analyst/FreeFormPointSet.java | 3 +-- src/main/java/com/conveyal/r5/analyst/Grid.java | 12 +++++------- 3 files changed, 7 insertions(+), 12 deletions(-) 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/FreeFormPointSet.java b/src/main/java/com/conveyal/r5/analyst/FreeFormPointSet.java index 4e9527be5..d0deaed13 100644 --- a/src/main/java/com/conveyal/r5/analyst/FreeFormPointSet.java +++ b/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. */ diff --git a/src/main/java/com/conveyal/r5/analyst/Grid.java b/src/main/java/com/conveyal/r5/analyst/Grid.java index 628ec31ec..dc6559330 100644 --- a/src/main/java/com/conveyal/r5/analyst/Grid.java +++ b/src/main/java/com/conveyal/r5/analyst/Grid.java @@ -122,7 +122,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 { @@ -615,7 +615,7 @@ public static List 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) { @@ -687,7 +687,7 @@ public static List 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 numericAttributes = reader.numericAttributes(); Set uniqueNumericAttributes = new HashSet<>(numericAttributes); if (uniqueNumericAttributes.size() != numericAttributes.size()) { @@ -764,14 +764,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++) { From 1d70428859a0393f7a646be541a65eb8d00aa8bc Mon Sep 17 00:00:00 2001 From: Andrew Byrd Date: Tue, 17 Oct 2023 00:18:13 +0200 Subject: [PATCH 3/6] Combine WGS84 to web Mercator functions WebMercatorGridPointSet methods now call static Grid functions with zoom. All functions are now using standard Java Math (not FastMath). --- .../java/com/conveyal/r5/analyst/Grid.java | 40 +++++++++++-------- .../r5/analyst/WebMercatorExtents.java | 7 +++- .../r5/analyst/WebMercatorGridPointSet.java | 37 ++++------------- .../r5/analyst/WebMercatorExtentsTest.java | 2 - 4 files changed, 38 insertions(+), 48 deletions(-) diff --git a/src/main/java/com/conveyal/r5/analyst/Grid.java b/src/main/java/com/conveyal/r5/analyst/Grid.java index dc6559330..ad9faf1ed 100644 --- a/src/main/java/com/conveyal/r5/analyst/Grid.java +++ b/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. @@ -462,6 +464,7 @@ public int featureCount() { /** 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; } @@ -474,12 +477,13 @@ public static int lonToPixel (double lon, int 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 x pixel number). Noninteger x pixel coordinates will return - * WGS84 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; } /** @@ -492,7 +496,7 @@ public static double pixelToCenterLon (int xPixel, 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); + final double latRad = toRadians(lat); return (1 - log(tan(latRad) + 1 / cos(latRad)) / Math.PI) * Math.pow(2, zoom - 1) * 256; } @@ -502,12 +506,16 @@ public static int latToPixel (double lat, int 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. + * 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 pixelToLat (double yPixel, int zoom) { - return FastMath.toDegrees(atan(sinh(Math.PI - (yPixel / 256d) / Math.pow(2, zoom) * 2 * Math.PI))); + final double tile = yPixel / 256d; + return toDegrees(atan(sinh(PI - tile * PI * 2 / pow(2, zoom)))); } /** diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java index b17a181c0..21a833d1b 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java +++ b/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java @@ -10,7 +10,12 @@ import java.util.Arrays; -import static com.conveyal.r5.analyst.Grid.*; +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; diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java index 46607a470..f8665421f 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java +++ b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java @@ -150,54 +150,33 @@ 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 - * of the world (not relative to this grid or to any particular tile). - * TODO This method could be reusable and static if zoom level was a parameter. And moved to WebMercatorExtents. + * See com.conveyal.r5.analyst.Grid#lonToPixel(double, int) */ public int lonToPixel (double lon) { - // factor of 256 is to get a pixel value not a tile number - return (int) ((lon + 180) / 360 * Math.pow(2, zoom) * 256); + return Grid.lonToPixel(lon, zoom); } /** - * Return the y pixel number containing the given latitude at this grid's zoom level, relative to the top edge of - * the world (not relative to this grid or to any particular tile). - * TODO This method could be reusable and static if zoom level was a parameter. And moved to WebMercatorExtents. + * See com.conveyal.r5.analyst.Grid#latToPixel(double, int) */ public int latToPixel (double lat) { - double invCos = 1 / Math.cos(Math.toRadians(lat)); - double tan = Math.tan(Math.toRadians(lat)); - double ln = Math.log(tan + invCos); - return (int) ((1 - ln / Math.PI) * Math.pow(2, zoom - 1) * 256); + return Grid.latToPixel(lat, zoom); } /** - * Return the longitude in degrees of the west edge of any pixel at the specified x coordinate relative to the left - * edge of the world (not relative to this grid or to any particular tile), at this grid's zoom level. - * TODO This method could be reusable and static if zoom level was a parameter. - * It should probably also be renamed to clarify that it doesn't return the center of the pixel. - * Another equivalent method is found in Grid, and should probably be merged with this one and WebMercatorExtents. + * See com.conveyal.r5.analyst.Grid#pixelToLon(double, int) */ public double pixelToLon (double x) { - return x / (Math.pow(2, zoom) * 256) * 360 - 180; + return Grid.pixelToLon(x, zoom); } /** - * Return the latitude in degrees of the north edge of any pixel at the specified y coordinate relative to the top - * edge of the world (not relative to this grid or to any particular tile), at this grid's zoom level. - * TODO This method could be reusable and static if zoom level was a parameter. - * It should probably also be renamed to clarify that it doesn't return the center of the pixel. - * Another equivalent method is found in Grid, and should probably be merged with this one and WebMercatorExtents. + * See com.conveyal.r5.analyst.Grid#pixelToLat(double, int) */ public double pixelToLat (double y) { - double tile = y / 256d; - return Math.toDegrees(Math.atan(Math.sinh(Math.PI - tile * Math.PI * 2 / Math.pow(2, zoom)))); + return Grid.pixelToLat(y, zoom); } @Override diff --git a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java index 5d0744556..2b4556154 100644 --- a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java +++ b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java @@ -4,8 +4,6 @@ import org.junit.jupiter.api.Test; import org.locationtech.jts.geom.Envelope; -import static org.junit.jupiter.api.Assertions.*; - class WebMercatorExtentsTest { /** From 6b05bae15c74fb973b1fcd967c4e640a1b2bce15 Mon Sep 17 00:00:00 2001 From: Andrew Byrd Date: Tue, 17 Oct 2023 01:38:15 +0200 Subject: [PATCH 4/6] Parameterize tests with zoom levels Also check stablity of single-pixel extents. --- .../r5/analyst/WebMercatorExtentsTest.java | 24 ++++++++++++------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java index 2b4556154..310166529 100644 --- a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java +++ b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java @@ -2,6 +2,8 @@ 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 { @@ -23,19 +25,20 @@ class WebMercatorExtentsTest { * 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 () { + @ParameterizedTest + @ValueSource(ints = {9, 10, 11}) + void wgsMercatorStability (int zoom) { Envelope wgsEnvelope = new Envelope(); wgsEnvelope.expandToInclude(10.22222222, 45.111111); wgsEnvelope.expandBy(0.1); - WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, 9); + 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, 9); + WebMercatorExtents webMercatorExtents2 = WebMercatorExtents.forTrimmedWgsEnvelope(wgsEnvelope, zoom); Assertions.assertEquals(webMercatorExtents2, webMercatorExtents); webMercatorExtents = webMercatorExtents2; } @@ -45,19 +48,24 @@ void wgsMercatorStability () { * 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 () { + @ParameterizedTest + @ValueSource(ints = {9, 10, 11}) + void singleCellExtents (int zoom) { Envelope wgsEnvelope = new Envelope(); wgsEnvelope.expandToInclude(10.22222222, 10.32222222); - WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, 9); + WebMercatorExtents webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, zoom); Assertions.assertEquals(1, webMercatorExtents.width); Assertions.assertEquals(1, webMercatorExtents.height); wgsEnvelope.expandBy(0.00001); - webMercatorExtents = WebMercatorExtents.forWgsEnvelope(wgsEnvelope, 9); + 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 From 3e9a86468a8b833eb1045860e26753b1231a6f50 Mon Sep 17 00:00:00 2001 From: Andrew Byrd Date: Tue, 17 Oct 2023 01:57:11 +0200 Subject: [PATCH 5/6] Clarify Javadoc --- .../r5/analyst/WebMercatorExtents.java | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java index 21a833d1b..d332d32f6 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java +++ b/src/main/java/com/conveyal/r5/analyst/WebMercatorExtents.java @@ -154,23 +154,26 @@ public static WebMercatorExtents forWgsEnvelope (Envelope wgsEnvelope, int zoom) 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 + * 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). - * 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. + * 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 to contain sets of - * points, it is probably better to deal with numerical imprecision by expanding the envelope slightly instead of - * contracting it. + * 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."); From d00317b9ed3f78361c2905401e2938efe1c248b4 Mon Sep 17 00:00:00 2001 From: Andrew Byrd Date: Tue, 17 Oct 2023 14:07:52 +0200 Subject: [PATCH 6/6] Add highest zoom level to tests, add comment --- .../java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java | 1 + .../java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java index f8665421f..f34023cb2 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java +++ b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java @@ -125,6 +125,7 @@ public double getLon(int i) { public TIntList getPointsInEnvelope (Envelope envelopeFixedDegrees) { // Convert fixed-degree envelope to floating, then to world-scale web Mercator pixels at this grid's zoom level. // This is not very DRY since we do something very similar in the constructor and elsewhere. + // These are the integer pixel numbers containing the envelope, so iteration below must be inclusive of the max. int west = lonToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMinX())); int east = lonToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMaxX())); int north = latToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMaxY())); diff --git a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java index 310166529..e0600cb6e 100644 --- a/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java +++ b/src/test/java/com/conveyal/r5/analyst/WebMercatorExtentsTest.java @@ -26,7 +26,7 @@ class WebMercatorExtentsTest { * digits (as when they are stored in a database or serialized over the wire). */ @ParameterizedTest - @ValueSource(ints = {9, 10, 11}) + @ValueSource(ints = {9, 10, 11, 12}) void wgsMercatorStability (int zoom) { Envelope wgsEnvelope = new Envelope(); wgsEnvelope.expandToInclude(10.22222222, 45.111111); @@ -49,7 +49,7 @@ void wgsMercatorStability (int zoom) { * 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}) + @ValueSource(ints = {9, 10, 11, 12}) void singleCellExtents (int zoom) { Envelope wgsEnvelope = new Envelope();