diff --git a/src/main/java/com/conveyal/r5/analyst/Grid.java b/src/main/java/com/conveyal/r5/analyst/Grid.java index d3d113ff1..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,15 +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 org.apache.commons.math3.util.FastMath.toRadians; +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. @@ -482,7 +483,7 @@ public static int lonToPixel (double lon, int zoom) { * 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; } /** @@ -505,20 +506,24 @@ public static int latToPixel (double lat, int 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); } /** diff --git a/src/main/java/com/conveyal/r5/analyst/IsochroneFeature.java b/src/main/java/com/conveyal/r5/analyst/IsochroneFeature.java deleted file mode 100644 index 893f20624..000000000 --- a/src/main/java/com/conveyal/r5/analyst/IsochroneFeature.java +++ /dev/null @@ -1,392 +0,0 @@ -package com.conveyal.r5.analyst; - -import com.conveyal.r5.common.GeometryUtils; -import com.google.common.collect.HashMultimap; -import com.google.common.collect.Multimap; -import org.locationtech.jts.geom.Coordinate; -import org.locationtech.jts.geom.LinearRing; -import org.locationtech.jts.geom.MultiPolygon; -import org.locationtech.jts.geom.Polygon; -import org.locationtech.jts.simplify.TopologyPreservingSimplifier; -import org.slf4j.Logger; -import org.slf4j.LoggerFactory; - -import java.io.File; -import java.io.FileWriter; -import java.io.Serializable; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Collection; -import java.util.Comparator; -import java.util.List; -import java.util.Map; -import java.util.stream.Collectors; - -/** - * This is similar to the IsochroneData class in OTP, and in fact for compatibility can be serialized to JSON and - * deserialized as such. However it uses a completely different algorithm. - * - * Although we have another separate implementaion of the marching squares contour line algorithm in Javascript, - * we're holding on to this one in case we need an R5 server to generate vector isochrones itself. - */ -public class IsochroneFeature implements Serializable { - private static final Logger LOG = LoggerFactory.getLogger(IsochroneFeature.class); - - private static final long serialVersionUID = 1L; - - // maximum number of vertices in a ring - public static final int MAX_RING_SIZE = 25000; - - /** The minimum ring size (to get rid of small rings). Should be at least 4 to ensure all rings are valid */ - public static final int MIN_RING_SIZE = 12; - - public MultiPolygon geometry; - public int cutoffSec; - - public IsochroneFeature () { /* deserialization */ } - - /** - * Create an isochrone for the given cutoff, using a Marching Squares algorithm. - * https://en.wikipedia.org/wiki/Marching_squares - */ - public IsochroneFeature (int cutoffSec, WebMercatorGridPointSet points, int[] times) { - // slightly hacky, but simple: set all of the times around the edges of the pointset to MAX_VALUE so that - // the isochrone never runs off the edge of the display. - // first, protective copy - times = Arrays.copyOf(times, times.length); - - for (int x = 0; x < points.width; x++) { - times[x] = Integer.MAX_VALUE; - times[(int) ((points.height - 1) * points.width) + x] = Integer.MAX_VALUE; - } - - for (int y = 0; y < points.height; y++) { - times[(int) points.width * y] = Integer.MAX_VALUE; - times[(int) points.width * (y + 1) - 1] = Integer.MAX_VALUE; - } - - LOG.debug("Making isochrone for {}sec", cutoffSec); - this.cutoffSec = cutoffSec; - // make contouring grid - byte[][] contour = new byte[(int) points.width - 1][(int) points.height - 1]; - for (int y = 0; y < points.height - 1; y++) { - for (int x = 0; x < points.width - 1; x++) { - boolean topLeft = times[(int) (points.width * y + x)] < cutoffSec; - boolean topRight = times[(int) (points.width * y + x + 1)] < cutoffSec; - boolean botLeft = times[(int) (points.width * (y + 1) + x)] < cutoffSec; - boolean botRight = times[(int) (points.width * (y + 1) + x + 1)] < cutoffSec; - - byte idx = 0; - - // TODO saddle points. Do we care? - - if (topLeft) idx |= 1 << 3; - if (topRight) idx |= 1 << 2; - if (botRight) idx |= 1 << 1; - if (botLeft) idx |= 1; - - contour[x][y] = idx; - } - } - - // create a geometry. For now not doing linear interpolation. Find a cell a line crosses through and - // follow that line. - List outerRings = new ArrayList<>(); - List innerRings = new ArrayList<>(); - boolean[][] found = new boolean[(int) points.width - 1][(int) points.height - 1]; - - for (int origy = 0; origy < points.height - 1; origy++) { - for (int origx = 0; origx < points.width - 1; origx++) { - int x = origx; - int y = origy; - - if (found[x][y]) continue; - - byte idx = contour[x][y]; - - // can't start at a saddle we don't know which way it goes - if (idx == 0 || idx == 5 || idx == 10 || idx == 15) continue; - - byte prevIdx = -1; - - List ring = new ArrayList<>(); - - // keep track of clockwise/counterclockwise orientation, see http://stackoverflow.com/questions/1165647 - int direction = 0; - // skip empty cells - int prevy = 0, prevx = 0; - Coordinate prevCoord = null; - CELLS: - while (true) { - idx = contour[x][y]; - - // check for intersecting rings, but know that saddles are supposed to self-intersect. - if (found[x][y] && idx != 5 && idx != 10) { - LOG.error("Ring crosses another ring (possibly itself). This cell has index {}, the previous cell has index {}.", idx, prevIdx); - break CELLS; - } - - found[x][y] = true; - - // follow line, keeping unfilled area to the left, which determines a direction - // this also means that we'll be able to figure out if something is a hole by - // the winding direction. - if (ring.size() >= MAX_RING_SIZE) { - LOG.error("Ring is too large, bailing"); - break CELLS; - } - - // save x values here, the next iteration may need to know what they were before we messed with them - // NB no bounds checking is performed below, but the next iteration of the loop will try to access contour[x][y] which - // will serve as a bounds check. - int startx = x; - int starty = y; - switch (idx) { - case 0: - LOG.error("Ran off outside of ring"); - break CELLS; - case 1: - x--; - break; - // NB: +y is down - case 2: - y++; - break; - case 3: - x--; - break; - case 4: - x++; - break; - case 5: - if (prevy > y) - // came from bottom - x++; - else if (prevy < y) - // came from top - x--; - else - LOG.error("Entered case 5 saddle point from wrong direction!"); - break; - case 6: - y++; - break; - case 7: - x--; - break; - case 8: - y--; - break; - case 9: - y--; - break; - case 10: - if (prevx < x) - // came from left - y++; - else if (prevx > x) - // came from right - y--; - else { - LOG.error("Entered case 10 saddle point from wrong direction."); - } - break; - case 11: - y--; - break; - case 12: - x++; - break; - case 13: - x++; - break; - case 14: - y++; - break; - case 15: - LOG.error("Ran off inside of ring"); - break CELLS; - } - - // figure out from whence we came - int topLeftTime = times[(int) points.width * y + x]; - int botLeftTime = times[(int) points.width * (y + 1) + x]; - int topRightTime = times[(int) points.width * y + x + 1]; - int botRightTime = times[(int) points.width * (y + 1) + x + 1]; - - double lat, lon; - - if (startx < x) { - // came from left - // will always be positive, if numerator is negative denominator will be as well. - double frac = (cutoffSec - topLeftTime) / (double) (botLeftTime - topLeftTime); - lat = points.pixelToLat(points.north + y + frac); - lon = points.pixelToLon(points.west + x); - } - else if (startx > x) { - // came from right - double frac = (cutoffSec - topRightTime) / (double) (botRightTime - topRightTime); - lat = points.pixelToLat(points.north + y + frac); - lon = points.pixelToLon(points.west + x + 1); - } - else if (starty < y) { - // came from top - double frac = (cutoffSec - topLeftTime) / (double) (topRightTime - topLeftTime); - lat = points.pixelToLat(points.north + y); - lon = points.pixelToLon(points.west + x + frac); - } - else { - // came from bottom - double frac = (cutoffSec - botLeftTime) / (botRightTime - botLeftTime); - lat = points.pixelToLat(points.north + y + 1); - lon = points.pixelToLon(points.west + x + frac); - } - - // keep track of winding direction - // http://stackoverflow.com/questions/1165647 - direction += (x - startx) * (y + starty); - - prevCoord = new Coordinate(lon, lat); - ring.add(prevCoord); - - // this shouldn't happen - if (x == startx && y == starty) { - LOG.error("Ring position did not update"); - break CELLS; - } - - // pass previous values to next iteration - prevIdx = idx; - prevx = startx; - prevy = starty; - - if (x == origx && y == origy) { - Coordinate end = ring.get(0); - ring.add(end); - - if (ring.size() > MIN_RING_SIZE) { - LinearRing lr = GeometryUtils.geometryFactory.createLinearRing(ring.toArray(new Coordinate[ring.size()])); - // direction less than 0 means clockwise (NB the y-axis is backwards), since value is to left it is an outer ring - if (direction > 0) { - // simplify so point in polygon test is tractable - lr = (LinearRing) TopologyPreservingSimplifier.simplify(lr, 1e-3); - outerRings.add(lr); - } else { - innerRings.add(lr); - } - } - - break CELLS; - } - } - } - } - - LOG.debug("{} components", outerRings.size()); - - Multimap holesForRing = HashMultimap.create(); - - // create polygons so we can test containment - Map polygonsForOuterRing = outerRings.stream().collect(Collectors.toMap( - r -> r, - r -> GeometryUtils.geometryFactory.createPolygon(r) - )); - - Map polygonsForInnerRing = innerRings.stream().collect(Collectors.toMap( - r -> r, - r -> GeometryUtils.geometryFactory.createPolygon(r) - )); - - // put the biggest ring first because most holes are in the biggest ring, reduces number of point in polygon tests below - outerRings.sort(Comparator.comparing(ring -> polygonsForOuterRing.get(ring).getArea()).reversed()); - - // get rid of tiny shells - - - LOG.info("Found {} outer rings and {} inner rings for cutoff {}m", polygonsForOuterRing.size(), polygonsForInnerRing.size(), cutoffSec / 60); - - int holeIdx = -1; - HOLES: for (Map.Entry hole : polygonsForInnerRing.entrySet()) { - holeIdx++; - - // get rid of tiny holes - if (hole.getValue().getArea() < 1e-6) continue; - - for (LinearRing ring : outerRings) { - // fine to test membership of first coordinate only since shells and holes are disjoint, and holes - // nest completely in shells - if (polygonsForOuterRing.get(ring).contains(hole.getKey().getPointN(0))) { - holesForRing.put(ring, hole.getKey()); - continue HOLES; - } - } - - LOG.warn("Found no fitting shell for isochrone hole {} at cutoff {}, dropping this hole.", holeIdx, cutoffSec); - } - - Polygon[] polygons = outerRings.stream().map(shell -> { - Collection holes = holesForRing.get(shell); - return GeometryUtils.geometryFactory.createPolygon(shell, holes.toArray(new LinearRing[holes.size()])); - }).toArray(s -> new Polygon[s]); - - // first geometry has to be an outer ring, but there may be multiple outer rings - this.geometry = GeometryUtils.geometryFactory.createMultiPolygon(polygons); - - LOG.debug("Done."); - } - - /** - * Debug code to draw an ascii-art isochrone. - */ - public void dumpContourGrid (byte[][] contour, String out) throws Exception { - FileWriter sb = new FileWriter(new File(out)); - for (int y = 0; y < contour[0].length; y++) { - for (int x = 0; x < contour.length; x++) { - switch(contour[x][y]) { - case 0: - sb.append(" "); - break; - case 15: - sb.append("XX"); - break; - case 1: - case 14: - sb.append("\\ "); - break; - case 2: - case 13: - sb.append(" /"); - break; - case 3: - case 12: - sb.append("--"); - break; - case 4: - case 11: - sb.append(" \\"); - break; - case 5: - sb.append("//"); - break; - case 10: - sb.append("\\\\"); - break; - case 6: - case 9: - sb.append("| "); - break; - case 7: - case 8: - sb.append("/ "); - break; - default: - sb.append("**"); - - } - } - sb.append("\n"); - } - - sb.close(); - } -} diff --git a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java index b6b614df7..bc5475c2c 100644 --- a/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java +++ b/src/main/java/com/conveyal/r5/analyst/WebMercatorGridPointSet.java @@ -10,6 +10,10 @@ import java.io.Serializable; +import static com.conveyal.r5.analyst.Grid.latToPixel; +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.conveyal.r5.common.GeometryUtils.checkWgsEnvelopeSize; import static com.conveyal.r5.streets.VertexStore.fixedDegreesToFloating; import static com.google.common.base.Preconditions.checkArgument; @@ -82,10 +86,10 @@ public WebMercatorGridPointSet (Envelope wgsEnvelope) { LOG.info("Creating WebMercatorGridPointSet with WGS84 extents {}", wgsEnvelope); checkWgsEnvelopeSize(wgsEnvelope, "grid point set"); this.zoom = DEFAULT_ZOOM; - int west = lonToPixel(wgsEnvelope.getMinX()); - int east = lonToPixel(wgsEnvelope.getMaxX()); - int north = latToPixel(wgsEnvelope.getMaxY()); - int south = latToPixel(wgsEnvelope.getMinY()); + int west = lonToPixel(wgsEnvelope.getMinX(), zoom); + int east = lonToPixel(wgsEnvelope.getMaxX(), zoom); + int north = latToPixel(wgsEnvelope.getMaxY(), zoom); + int south = latToPixel(wgsEnvelope.getMinY(), zoom); this.west = west; this.north = north; this.height = south - north; @@ -112,23 +116,23 @@ public double sumTotalOpportunities () { @Override public double getLat(int i) { long y = i / this.width + this.north; - return pixelToLat(y); + return pixelToLat(y, zoom); } @Override public double getLon(int i) { long x = i % this.width + this.west; - return pixelToLon(x); + return pixelToLon(x, zoom); } @Override 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. - int west = lonToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMinX())); - int east = lonToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMaxX())); - int north = latToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMaxY())); - int south = latToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMinY())); + int west = lonToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMinX()), zoom); + int east = lonToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMaxX()), zoom); + int north = latToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMaxY()), zoom); + int south = latToPixel(fixedDegreesToFloating(envelopeFixedDegrees.getMinY()), zoom); // Make the envelope's pixel values relative to the edges of this WebMercatorGridPointSet, rather than // absolute world-scale coordinates at this zoom level. west -= this.west; @@ -149,53 +153,6 @@ public TIntList getPointsInEnvelope (Envelope envelopeFixedDegrees) { return pointsInEnvelope; } - // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Mathematics - - /** - * 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. - */ - 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 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. - */ - 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 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. - */ - public double pixelToLon (double x) { - return x / (Math.pow(2, zoom) * 256) * 360 - 180; - } - - /** - * 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. - */ - 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)))); - } - @Override public double getOpportunityCount (int i) { // FIXME just counting the points for now, return counts @@ -215,8 +172,8 @@ public int getPointIndexContaining (Coordinate latLon) { } public int getPointIndexContaining (double lon, double lat) { - int x = lonToPixel(lon) - west; - int y = latToPixel(lat) - north; + int x = lonToPixel(lon, zoom) - west; + int y = latToPixel(lat, zoom) - north; return y * width + x; } diff --git a/src/test/java/com/conveyal/r5/analyst/GridTest.java b/src/test/java/com/conveyal/r5/analyst/GridTest.java index 03149d1c2..7e5bf62a4 100644 --- a/src/test/java/com/conveyal/r5/analyst/GridTest.java +++ b/src/test/java/com/conveyal/r5/analyst/GridTest.java @@ -177,5 +177,17 @@ private static void assertArrayEquals ( assertEquals(aSum, bSum, gridSumTolerance); } } + + /** Test that latitude/longitude to pixel conversions are correct */ + @Test + public void testLatLonPixelConversions () { + final int ZOOM = WebMercatorGridPointSet.DEFAULT_ZOOM; + for (double lat : new double [] { -75, -25, 0, 25, 75 }) { + assertEquals(lat, Grid.pixelToLat(Grid.latToPixel(lat, ZOOM), ZOOM), 1e-2); + } + for (double lon : new double [] { -175, -90, 0, 90, 175}) { + assertEquals(lon, Grid.pixelToLon(Grid.lonToPixel(lon, ZOOM), ZOOM), 1e-2); + } + } } diff --git a/src/test/java/com/conveyal/r5/analyst/WebMercatorGridPointSetTest.java b/src/test/java/com/conveyal/r5/analyst/WebMercatorGridPointSetTest.java index 5f85ae2b7..a3bd481c7 100644 --- a/src/test/java/com/conveyal/r5/analyst/WebMercatorGridPointSetTest.java +++ b/src/test/java/com/conveyal/r5/analyst/WebMercatorGridPointSetTest.java @@ -20,23 +20,6 @@ */ public class WebMercatorGridPointSetTest { - /** Test that latitude/longitude to pixel conversions are correct */ - @Test - public void testLatLonPixelConversions () { - // offsets and width/height not needed for calculations - WebMercatorGridPointSet ps = new WebMercatorGridPointSet( - WebMercatorGridPointSet.DEFAULT_ZOOM, 0,0, 100, 100, null - ); - - for (double lat : new double [] { -75, -25, 0, 25, 75 }) { - assertEquals(lat, ps.pixelToLat(ps.latToPixel(lat)), 1e-2); - } - - for (double lon : new double [] { -175, -90, 0, 90, 175}) { - assertEquals(lon, ps.pixelToLon(ps.lonToPixel(lon)), 1e-2); - } - } - /** Test that we can find pixel numbers for envelopes intersecting a gridded PointSet. */ @Test public void testPixelsInEnvelope () {