diff --git a/docs/Configuration.md b/docs/Configuration.md index 4f2fd2fc47e..a211edc3ec6 100644 --- a/docs/Configuration.md +++ b/docs/Configuration.md @@ -84,7 +84,8 @@ config key | description | value type | value default | notes `matchBusRoutesToStreets` | Based on GTFS shape data, guess which OSM streets each bus runs on to improve stop linking | boolean | false | `fetchElevationUS` | Download US NED elevation data and apply it to the graph | boolean | false | `elevationBucket` | If specified, download NED elevation tiles from the given AWS S3 bucket | object | null | provide an object with `accessKey`, `secretKey`, and `bucketName` for AWS S3 -`elevationUnitMultiplier` | Specify a multiplier to convert elevation units from source to meters | double | 1.0 | see [Elevation unit conversion](#elevation-unit-conversion) +`readCachedElevations` | If true, reads in pre-calculated elevation data. | boolean | true | see [Elevation Data Calculation Optimizations](#elevation-data-calculation-optimizations) +`writeCachedElevations` | If true, writes the calculated elevation data. | boolean | false | see [Elevation Data Calculation Optimizations](#elevation-data-calculation-optimizations) `fares` | A specific fares service to use | object | null | see [fares configuration](#fares-configuration) `osmNaming` | A custom OSM namer to use | object | null | see [custom naming](#custom-naming) `osmWayPropertySet` | Custom OSM way properties | string | `default` | options: `default`, `norway` @@ -197,6 +198,31 @@ You can configure it as follows in `build-config.json`: } ``` +### Geoid Difference + +With some elevation data, the elevation values are specified as relative to the a geoid (irregular estimate of mean sea level). See [issue #2301](https://github.com/opentripplanner/OpenTripPlanner/issues/2301) for detailed discussion of this. In these cases, it is necessary to also add this geoid value onto the elevation value to get the correct result. OTP can automatically calculate these values in one of two ways. + +The first way is to use the geoid difference value that is calculated once at the center of the graph. This value is returned in each trip plan response in the [ElevationMetadata](http://otp-docs.ibi-transit.com/api/json_ElevationMetadata.html) field. Using a single value can be sufficient for smaller OTP deployments, but might result in incorrect values at the edges of larger OTP deployments. If your OTP instance uses this, it is recommended to set a default request value in the `router-config.json` file as follows: + +```JSON +// router-config.json +{ + "routingDefaults": { + "geoidElevation ": true + } +} +``` + +The second way is to precompute these geoid difference values at a more granular level and include them when calculating elevations for each sampled point along each street edge. In order to speed up calculations, the geoid difference values are calculated and cached using only 2 significant digits of GPS coordinates. This is more than enough detail for most regions of the world and should result in less than one meter of difference in areas that have large changes in geoid difference values. To enable this, include the following in the `build-config.json` file: + +```JSON +// build-config.json +{ + "includeEllipsoidToGeoidDifference": true +} +``` + +If the geoid difference values are precomputed, be careful to not set the routing resource value of `geoidElevation` to true in order to avoid having the graph-wide geoid added again to all elevation values in the relevant street edges in responses. ### Other raster elevation data @@ -218,6 +244,28 @@ order as the above-mentioned SRTM data, which is also the default for the popula DEM files(USGS DEM) is not supported by OTP, but can be converted to GeoTIFF with tools like [GDAL](http://www.gdal.org/). Use `gdal_merge.py -o merged.tiff *.dem` to merge a set of `dem` files into one `tif` file. +### Elevation Data Calculation Optimizations + +Calculating elevations on all StreetEdges can take a dramatically long time. In a very large graph build for multiple Northeast US states, the time it took to download the elevation data and calculate all of the elevations took 5,509 seconds (roughly 1.5 hours). + +If you are using cloud computing for your OTP instances, it is recommended to create prebuilt images that contain the elevation data you need. This will save time because all of the data won't need to be downloaded. + +However, the bulk of the time will still be spent calculating elevations for all of the street edges. Therefore, a further optimization can be done to calculate and save the elevation data during a graph build and then save it for future use. + +In order to write out the precalculated elevation data, add this to your `build-config.json` file: + +```JSON +// build-config.json +{ + "writeCachedElevations": true +} +``` + +After building the graph, a file called `cached_elevations.obj` will be written to the cache directory. By default, this file is not written during graph builds. There is also a graph build parameter called `readCachedElevations` which is set to `true` by default. + +In graph builds, the elevation module will attempt to read the `cached_elevations.obj` file from the cache directory. The cache directory defaults to `/var/otp/cache`, but this can be overriden via the CLI argument `--cache `. For the same graph build for multiple Northeast US states, the time it took with using this predownloaded and precalculated data became 543.7 seconds (roughly 9 minutes). + +The cached data is a lookup table where the coordinate sequences of respective street edges are used as keys for calculated data. Therefore, it is expected that over time various edits to OpenStreetMap will cause this cached data to become stale and not include new OSM ways. Therefore, periodic update of this cached data is recommended. ## Fares configuration diff --git a/src/main/java/org/opentripplanner/graph_builder/GraphBuilder.java b/src/main/java/org/opentripplanner/graph_builder/GraphBuilder.java index c70b637afb7..6f1dcd5bc56 100644 --- a/src/main/java/org/opentripplanner/graph_builder/GraphBuilder.java +++ b/src/main/java/org/opentripplanner/graph_builder/GraphBuilder.java @@ -35,6 +35,7 @@ import java.io.File; import java.io.IOException; +import java.nio.file.Paths; import java.util.ArrayList; import java.util.HashMap; import java.util.List; @@ -179,6 +180,7 @@ public static GraphBuilder forDirectory(CommandLineParameters params, File dir) LOG.error("'{}' is not a readable directory.", dir); return null; } + // Find and parse config files first to reveal syntax errors early without waiting for graph build. builderConfig = OTPMain.loadJson(new File(dir, BUILDER_CONFIG_FILENAME)); GraphBuilderParameters builderParams = new GraphBuilderParameters(builderConfig); @@ -272,30 +274,34 @@ public static GraphBuilder forDirectory(CommandLineParameters params, File dir) graphBuilder.addModule(streetLinkerModule); // Load elevation data and apply it to the streets. // We want to do run this module after loading the OSM street network but before finding transfers. + ElevationGridCoverageFactory gcf = null; if (builderParams.elevationBucket != null) { // Download the elevation tiles from an Amazon S3 bucket S3BucketConfig bucketConfig = builderParams.elevationBucket; File cacheDirectory = new File(params.cacheDirectory, "ned"); DegreeGridNEDTileSource awsTileSource = new DegreeGridNEDTileSource(); - awsTileSource = new DegreeGridNEDTileSource(); awsTileSource.awsAccessKey = bucketConfig.accessKey; awsTileSource.awsSecretKey = bucketConfig.secretKey; awsTileSource.awsBucketName = bucketConfig.bucketName; - NEDGridCoverageFactoryImpl gcf = new NEDGridCoverageFactoryImpl(cacheDirectory); - gcf.tileSource = awsTileSource; - GraphBuilderModule elevationBuilder = new ElevationModule(gcf); - graphBuilder.addModule(elevationBuilder); + gcf = new NEDGridCoverageFactoryImpl(cacheDirectory, awsTileSource); } else if (builderParams.fetchElevationUS) { // Download the elevation tiles from the official web service File cacheDirectory = new File(params.cacheDirectory, "ned"); - ElevationGridCoverageFactory gcf = new NEDGridCoverageFactoryImpl(cacheDirectory); - GraphBuilderModule elevationBuilder = new ElevationModule(gcf); - graphBuilder.addModule(elevationBuilder); + gcf = new NEDGridCoverageFactoryImpl(cacheDirectory); } else if (demFile != null) { // Load the elevation from a file in the graph inputs directory - ElevationGridCoverageFactory gcf = new GeotiffGridCoverageFactoryImpl(demFile); - GraphBuilderModule elevationBuilder = new ElevationModule(gcf); - graphBuilder.addModule(elevationBuilder); + gcf = new GeotiffGridCoverageFactoryImpl(demFile); + } + if (gcf != null) { + graphBuilder.addModule( + new ElevationModule( + gcf, + params.cacheDirectory, + builderParams.readCachedElevations, + builderParams.writeCachedElevations, + builderParams.includeEllipsoidToGeoidDifference + ) + ); } if ( hasGTFS ) { // The stops can be linked to each other once they are already linked to the street network. diff --git a/src/main/java/org/opentripplanner/graph_builder/module/ned/DegreeGridNEDTileSource.java b/src/main/java/org/opentripplanner/graph_builder/module/ned/DegreeGridNEDTileSource.java index f3d887447b5..678d11c98a3 100644 --- a/src/main/java/org/opentripplanner/graph_builder/module/ned/DegreeGridNEDTileSource.java +++ b/src/main/java/org/opentripplanner/graph_builder/module/ned/DegreeGridNEDTileSource.java @@ -33,8 +33,6 @@ public class DegreeGridNEDTileSource implements NEDTileSource { private static Logger log = LoggerFactory.getLogger(DegreeGridNEDTileSource.class); - private Graph graph; - private File cacheDirectory; public String awsAccessKey; @@ -43,18 +41,15 @@ public class DegreeGridNEDTileSource implements NEDTileSource { public String awsBucketName; - @Override - public void setGraph(Graph graph) { - this.graph = graph; - } + private List nedTiles; @Override - public void setCacheDirectory(File cacheDirectory) { - this.cacheDirectory = cacheDirectory; + public List getNEDTiles() { + return nedTiles; } - @Override - public List getNEDTiles() { + @Override public void fetchData(Graph graph, File cacheDirectory) { + this.cacheDirectory = cacheDirectory; HashSet> tiles = new HashSet>(); @@ -67,9 +62,12 @@ public List getNEDTiles() { for (P2 tile : tiles) { int x = tile.first - 1; int y = tile.second + 1; - paths.add(getPathToTile(x, y)); + File tilePath = getPathToTile(x, y); + if (tilePath != null) { + paths.add(tilePath); + } } - return paths; + nedTiles = paths; } private String formatLatLon(int x, int y) { @@ -97,14 +95,15 @@ private File getPathToTile(int x, int y) { path.getParentFile().mkdirs(); if (awsAccessKey == null || awsSecretKey == null) { - throw new RuntimeException("Cannot download NED tiles from S3: awsAccessKey or awsSecretKey properties are not set"); + throw new RuntimeException( + "Cannot download NED tiles from S3: awsAccessKey or awsSecretKey properties are not set"); } log.info("Downloading NED degree tile " + path); // download the file from S3. AWSCredentials awsCredentials = new AWSCredentials(awsAccessKey, awsSecretKey); + String key = formatLatLon(x, y) + ".tiff"; try { S3Service s3Service = new RestS3Service(awsCredentials); - String key = formatLatLon(x, y) + ".tiff"; S3Object object = s3Service.getObject(awsBucketName, key); InputStream istream = object.getDataInputStream(); @@ -120,18 +119,10 @@ private File getPathToTile(int x, int y) { } ostream.close(); istream.close(); - } catch (S3ServiceException e) { - path.deleteOnExit(); - throw new RuntimeException(e); - } catch (ServiceException e) { - path.deleteOnExit(); - throw new RuntimeException(e); - } catch (FileNotFoundException e) { - path.deleteOnExit(); - throw new RuntimeException(e); - } catch (IOException e) { + } catch (ServiceException | IOException e) { + log.error("Error downloading tile {}! Error: {}.", key, e.getMessage()); path.deleteOnExit(); - throw new RuntimeException(e); + return null; } return path; } diff --git a/src/main/java/org/opentripplanner/graph_builder/module/ned/ElevationModule.java b/src/main/java/org/opentripplanner/graph_builder/module/ned/ElevationModule.java index 0bacd8536e8..bccbdaf97b2 100644 --- a/src/main/java/org/opentripplanner/graph_builder/module/ned/ElevationModule.java +++ b/src/main/java/org/opentripplanner/graph_builder/module/ned/ElevationModule.java @@ -2,10 +2,10 @@ import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Geometry; -import org.geotools.coverage.grid.GridCoverage2D; -import org.geotools.coverage.grid.Interpolator2D; import org.geotools.geometry.DirectPosition2D; import org.opengis.coverage.Coverage; +import org.opengis.coverage.PointOutsideCoverageException; +import org.opengis.referencing.operation.TransformException; import org.opentripplanner.common.geometry.GeometryUtils; import org.opentripplanner.common.geometry.PackedCoordinateSequence; import org.opentripplanner.common.geometry.SphericalDistanceLibrary; @@ -14,7 +14,6 @@ import org.opentripplanner.graph_builder.annotation.Graphwide; import org.opentripplanner.graph_builder.module.GraphBuilderModuleSummary; import org.opentripplanner.graph_builder.module.GraphBuilderTaskSummary; -import org.opentripplanner.graph_builder.module.extra_elevation_data.ElevationPoint; import org.opentripplanner.graph_builder.services.GraphBuilderModule; import org.opentripplanner.graph_builder.services.ned.ElevationGridCoverageFactory; import org.opentripplanner.routing.edgetype.StreetEdge; @@ -22,38 +21,62 @@ import org.opentripplanner.routing.graph.Edge; import org.opentripplanner.routing.graph.Graph; import org.opentripplanner.routing.graph.Vertex; +import org.opentripplanner.util.PolylineEncoder; import org.slf4j.Logger; import org.slf4j.LoggerFactory; -import javax.media.jai.InterpolationBilinear; -import java.util.ArrayList; +import java.io.BufferedOutputStream; +import java.io.File; +import java.io.FileInputStream; +import java.io.FileOutputStream; +import java.io.IOException; +import java.io.ObjectInputStream; +import java.io.ObjectOutputStream; +import java.nio.file.Files; import java.util.Arrays; import java.util.HashMap; import java.util.HashSet; import java.util.LinkedList; import java.util.List; +import java.util.concurrent.ConcurrentHashMap; +import java.util.concurrent.ForkJoinPool; +import java.util.concurrent.TimeUnit; +import java.util.concurrent.atomic.AtomicInteger; + +import static org.opentripplanner.util.ElevationUtils.computeEllipsoidToGeoidDifference; /** - * {@link org.opentripplanner.graph_builder.services.GraphBuilderModule} plugin that applies elevation data to street data that has already - * been loaded into a (@link Graph}, creating elevation profiles for each Street encountered - * in the Graph. Depending on the {@link ElevationGridCoverageFactory} specified - * this could be auto-downloaded and cached National Elevation Dataset (NED) raster data or - * a GeoTIFF file. The elevation profiles are stored as {@link PackedCoordinateSequence} objects, - * where each (x,y) pair represents one sample, with the x-coord representing the distance along - * the edge measured from the start, and the y-coord representing the sampled elevation at that - * point (both in meters). + * {@link org.opentripplanner.graph_builder.services.GraphBuilderModule} plugin that applies elevation data to street + * data that has already been loaded into a (@link Graph}, creating elevation profiles for each Street encountered + * in the Graph. Data sources that could be used include auto-downloaded and cached National Elevation Dataset (NED) + * raster data or a GeoTIFF file. The elevation profiles are stored as {@link PackedCoordinateSequence} objects, where + * each (x,y) pair represents one sample, with the x-coord representing the distance along the edge measured from the + * start, and the y-coord representing the sampled elevation at that point (both in meters). */ public class ElevationModule implements GraphBuilderModule { private static final Logger log = LoggerFactory.getLogger(ElevationModule.class); - private ElevationGridCoverageFactory gridCoverageFactory; - - private Coverage coverage; - - // Keep track of the proportion of elevation fetch operations that fail so we can issue warnings. - private int nPointsEvaluated = 0; - private int nPointsOutsideDEM = 0; + /** The elevation data to be used in calculating elevations. */ + private final ElevationGridCoverageFactory gridCoverageFactory; + /* Whether or not to attempt reading in a file of cached elevations */ + private final boolean readCachedElevations; + /* Whether or not to attempt writing out a file of cached elevations */ + private final boolean writeCachedElevations; + /* The file of cached elevations */ + private final File cachedElevationsFile; + /* Whether or not to include geoid difference values in individual elevation calculations */ + private final boolean includeEllipsoidToGeoidDifference; + + private HashMap cachedElevations; + + // Keep track of the proportion of elevation fetch operations that fail so we can issue warnings. AtomicInteger is + // used to provide thread-safe updating capabilities. + private AtomicInteger nEdgesProcessed = new AtomicInteger(0); + private final AtomicInteger nPointsEvaluated = new AtomicInteger(0); + private final AtomicInteger nPointsOutsideDEM = new AtomicInteger(0); + /** keeps track of the total amount of elevation edges for logging purposes */ + private int totalElevationEdges = Integer.MAX_VALUE; /** * The distance between samples in meters. Defaults to 10m, the approximate resolution of 1/3 @@ -61,10 +84,33 @@ public class ElevationModule implements GraphBuilderModule { */ private double distanceBetweenSamplesM = 10; - public ElevationModule() { /* This makes me a "bean" */ }; - + /** the graph being built */ + private Graph graph; + + /** A concurrent hashmap used for storing geoid difference values at various coordinates */ + private final ConcurrentHashMap geoidDifferenceCache = new ConcurrentHashMap<>(); + + // used only for testing purposes public ElevationModule(ElevationGridCoverageFactory factory) { - this.setGridCoverageFactory(factory); + gridCoverageFactory = factory; + cachedElevationsFile = null; + readCachedElevations = false; + writeCachedElevations = false; + includeEllipsoidToGeoidDifference = true; + } + + public ElevationModule( + ElevationGridCoverageFactory factory, + File cacheDirectory, + boolean readCachedElevations, + boolean writeCachedElevations, + boolean includeEllipsoidToGeoidDifference + ) { + gridCoverageFactory = factory; + cachedElevationsFile = new File(cacheDirectory, "cached_elevations.obj"); + this.readCachedElevations = readCachedElevations; + this.writeCachedElevations = writeCachedElevations; + this.includeEllipsoidToGeoidDifference = includeEllipsoidToGeoidDifference; } public List provides() { @@ -74,30 +120,28 @@ public List provides() { public List getPrerequisites() { return Arrays.asList("streets"); } - - public void setGridCoverageFactory(ElevationGridCoverageFactory factory) { - gridCoverageFactory = factory; - } - - public void setDistanceBetweenSamplesM(double distance) { - distanceBetweenSamplesM = distance; - } @Override public void buildGraph(Graph graph, GraphBuilderModuleSummary graphBuilderModuleSummary) { + this.graph = graph; GraphBuilderTaskSummary demPrepareTask = graphBuilderModuleSummary.addSubTask( "fetch and prepare elevation data" ); log.info(demPrepareTask.start()); - gridCoverageFactory.setGraph(graph); - Coverage gridCov = gridCoverageFactory.getGridCoverage(); - - // If gridCov is a GridCoverage2D, apply a bilinear interpolator. Otherwise, just use the - // coverage as is (note: UnifiedGridCoverages created by NEDGridCoverageFactoryImpl handle - // interpolation internally) - coverage = (gridCov instanceof GridCoverage2D) ? Interpolator2D.create( - (GridCoverage2D) gridCov, new InterpolationBilinear()) : gridCov; + gridCoverageFactory.fetchData(graph); + + // try to load in the cached elevation data + if (readCachedElevations) { + try { + ObjectInputStream in = new ObjectInputStream(new FileInputStream(cachedElevationsFile)); + cachedElevations = (HashMap) in.readObject(); + log.info("Cached elevation data loaded into memory!"); + } catch (IOException | ClassNotFoundException e) { + log.warn(graph.addBuilderAnnotation(new Graphwide( + String.format("Cached elevations file could not be read in due to error: %s!", e.getMessage())))); + } + } log.info(demPrepareTask.finish()); GraphBuilderTaskSummary setElevationsFromDEMTask = graphBuilderModuleSummary.addSubTask( @@ -105,46 +149,134 @@ public void buildGraph(Graph graph, GraphBuilderModuleSummary graphBuilderModule ); log.info(setElevationsFromDEMTask.start()); - List edgesWithElevation = new ArrayList(); - int nProcessed = 0; - int nTotal = graph.countEdges(); + // Multithread elevation calculations + ForkJoinPool forkJoinPool = new ForkJoinPool(); + + // For unknown reasons, the interpolation of heights at coordinates is a synchronized method in the commonly + // used Interpolator2D class. Therefore, it is critical to use a dedicated Coverage instance for each thread to + // avoid other threads waiting for a lock to be released on the Coverage instance. This concurrent HashMap will + // store these thread-specific Coverage instances. + ConcurrentHashMap coveragesForThread = new ConcurrentHashMap<>(); + + // At first, set the totalElevationEdges to the total number of edges in the graph. + totalElevationEdges = graph.countEdges(); + List streetsWithElevationEdges = new LinkedList<>(); for (Vertex gv : graph.getVertices()) { for (Edge ee : gv.getOutgoing()) { if (ee instanceof StreetWithElevationEdge) { - StreetWithElevationEdge edgeWithElevation = (StreetWithElevationEdge) ee; - processEdge(graph, edgeWithElevation); - if (edgeWithElevation.getElevationProfile() != null && !edgeWithElevation.isElevationFlattened()) { - edgesWithElevation.add(edgeWithElevation); - } - nProcessed += 1; - if (nProcessed % 50000 == 0) { - log.info("set elevation on {}/{} edges", nProcessed, nTotal); - double failurePercentage = nPointsOutsideDEM / nPointsEvaluated * 100; - if (failurePercentage > 50) { - log.warn(graph.addBuilderAnnotation(new Graphwide( - String.format( - "Fetching elevation failed at %d/%d points (%d%%)", - nPointsOutsideDEM, nPointsEvaluated, failurePercentage - ) - ))); - log.warn("Elevation is missing at a large number of points. DEM may be for the wrong region. " + - "If it is unprojected, perhaps the axes are not in (longitude, latitude) order."); - } - } + forkJoinPool.submit(new ProcessEdgeTask((StreetWithElevationEdge) ee, coveragesForThread)); + streetsWithElevationEdges.add((StreetWithElevationEdge) ee); } } } + // update this value to the now-known amount of edges that are StreetWithElevation edges + totalElevationEdges = streetsWithElevationEdges.size(); + + // shutdown the forkJoinPool and wait until all tasks are finished. If this takes longer than 1 day, give up. + forkJoinPool.shutdown(); + try { + forkJoinPool.awaitTermination(1, TimeUnit.DAYS); + } catch (InterruptedException e) { + log.warn(graph.addBuilderAnnotation(new Graphwide("Multi-threaded elevation calculations timed-out!"))); + } + + double failurePercentage = nPointsOutsideDEM.get() / nPointsEvaluated.get() * 100; + if (failurePercentage > 50) { + log.warn(graph.addBuilderAnnotation(new Graphwide( + String.format( + "Fetching elevation failed at %d/%d points (%d%%)", + nPointsOutsideDEM, nPointsEvaluated, failurePercentage + ) + ))); + log.warn("Elevation is missing at a large number of points. DEM may be for the wrong region. " + + "If it is unprojected, perhaps the axes are not in (longitude, latitude) order."); + } + // iterate again to find edges that had elevation calculated. This is done here instead of in the forkJoinPool + // to avoid thread locking for writes to a synchronized list + LinkedList edgesWithCalculatedElevations = new LinkedList<>(); + for (StreetWithElevationEdge edgeWithElevation : streetsWithElevationEdges) { + if (edgeWithElevation.hasPackedElevationProfile() && !edgeWithElevation.isElevationFlattened()) { + edgesWithCalculatedElevations.add(edgeWithElevation); + } + } + + if (writeCachedElevations) { + // write information from edgesWithElevation to a new cache file for subsequent graph builds + HashMap newCachedElevations = new HashMap<>(); + for (StreetEdge streetEdge : edgesWithCalculatedElevations) { + newCachedElevations.put(PolylineEncoder.createEncodings(streetEdge.getGeometry()).getPoints(), + streetEdge.getElevationProfile()); + } + try { + ObjectOutputStream out = new ObjectOutputStream( + new BufferedOutputStream(new FileOutputStream(cachedElevationsFile))); + out.writeObject(newCachedElevations); + out.close(); + } catch (IOException e) { + log.error(e.getMessage()); + log.error(graph.addBuilderAnnotation(new Graphwide("Failed to write cached elevation file!"))); + } + } log.info(setElevationsFromDEMTask.finish()); GraphBuilderTaskSummary missingElevationsTask = graphBuilderModuleSummary.addSubTask( "calculate missing elevations" ); log.info(missingElevationsTask.start()); - assignMissingElevations(graph, edgesWithElevation); + assignMissingElevations(graph, edgesWithCalculatedElevations); log.info(missingElevationsTask.finish()); } + /** + * A runnable that contains the relevant info for executing a process edge operation in a particular thread. + */ + private class ProcessEdgeTask implements Runnable { + private final StreetWithElevationEdge swee; + private final ConcurrentHashMap coveragesForThread; + + public ProcessEdgeTask(StreetWithElevationEdge swee, ConcurrentHashMap coveragesForThread) { + this.swee = swee; + this.coveragesForThread = coveragesForThread; + } + + @Override public void run() { + processEdge(swee, getThreadSpecificCoverageInstance()); + int curNumProcessed = nEdgesProcessed.addAndGet(1); + if (curNumProcessed % 50000 == 0) { + log.info("set elevation on {}/{} edges", curNumProcessed, totalElevationEdges); + } + } + + /** + * Get the thread-specific Coverage instance to avoid multiple threads waiting for a lock to be released on a + * Interpolator2D instance. + */ + private Coverage getThreadSpecificCoverageInstance () { + long currentThreadId = Thread.currentThread().getId(); + Coverage threadSpecificCoverage = coveragesForThread.get(currentThreadId); + if (threadSpecificCoverage == null) { + // Synchronize the creation of the Thread-specific Coverage instances to avoid potential locks that + // could arise from downstream classes that have synchronized methods. + synchronized (coveragesForThread) { + // Get a new Coverage instance from the module's ElevationGridCoverageFactory. + threadSpecificCoverage = gridCoverageFactory.getGridCoverage(); + // The Coverage instance relies on some synchronized static methods shared across all threads that + // can cause deadlocks if not fully initialized. Therefore, make a single request for the first + // point on the edge to initialize these other items. + Coordinate firstEdgeCoord = swee.getGeometry().getCoordinates()[0]; + double[] dummy = new double[1]; + threadSpecificCoverage.evaluate( + new DirectPosition2D(GeometryUtils.WGS84_XY, firstEdgeCoord.x, firstEdgeCoord.y), + dummy + ); + } + coveragesForThread.put(currentThreadId, threadSpecificCoverage); + } + return threadSpecificCoverage; + } + } + class ElevationRepairState { /* This uses an intuitionist approach to elevation inspection */ public StreetEdge backEdge; @@ -173,7 +305,7 @@ public ElevationRepairState(StreetEdge backEdge, ElevationRepairState backState, */ private void assignMissingElevations( Graph graph, - List edgesWithElevation + LinkedList edgesWithElevation ) { log.debug("Assigning missing elevations"); @@ -346,116 +478,123 @@ private void assignMissingElevations( /** * Processes a single street edge, creating and assigning the elevation profile. - * + * * @param ee the street edge - * @param graph the graph (used only for error handling) + * @param coverage the specific Coverage instance to use in order to avoid competition between threads */ - private void processEdge(Graph graph, StreetWithElevationEdge ee) { - if (ee.getElevationProfile() != null) { + private void processEdge(StreetWithElevationEdge ee, Coverage coverage) { + if (ee.hasPackedElevationProfile()) { return; /* already set up */ } Geometry g = ee.getGeometry(); - Coordinate[] coords = g.getCoordinates(); - List coordList = new LinkedList(); - - // calculate the total edge length in meters - double edgeLenM = 0; - for (int i = 0; i < coords.length - 1; i++) { - edgeLenM += SphericalDistanceLibrary.distance(coords[i].y, coords[i].x, coords[i + 1].y, - coords[i + 1].x); + // first try to find a cached value if possible + if (cachedElevations != null) { + PackedCoordinateSequence coordinateSequence = cachedElevations.get( + PolylineEncoder.createEncodings(g).getPoints() + ); + // found a cached value! + if (coordinateSequence != null) { + setEdgeElevationProfile(ee, coordinateSequence, graph); + return; + } } - // initial sample (x = 0) - coordList.add(new Coordinate(0, getElevation(coords[0]))); - - // loop for edge-internal samples - for (double x = distanceBetweenSamplesM; x < edgeLenM; x += distanceBetweenSamplesM) { - // avoid final-segment samples less than half the distance between samples: - if (edgeLenM - x < distanceBetweenSamplesM / 2) { - break; + // did not find a cached value, calculate + // If any of the coordinates throw an error when trying to lookup their value, immediately bail and do not + // process the elevation on the edge + try { + Coordinate[] coords = g.getCoordinates(); + + List coordList = new LinkedList(); + + // initial sample (x = 0) + coordList.add(new Coordinate(0, getElevation(coverage, coords[0]))); + + // iterate through coordinates calculating the edge length and creating intermediate elevation coordinates at + // the regularly specified interval + double edgeLenM = 0; + double sampleDistance = distanceBetweenSamplesM; + double previousDistance = 0; + double x1 = coords[0].x, y1 = coords[0].y, x2, y2; + for (int i = 0; i < coords.length - 1; i++) { + x2 = coords[i + 1].x; + y2 = coords[i + 1].y; + double curSegmentDistance = SphericalDistanceLibrary.distance(y1, x1, y2, x2); + edgeLenM += curSegmentDistance; + while (edgeLenM > sampleDistance) { + // if current edge length is longer than the current sample distance, insert new elevation coordinates + // as needed until sample distance has caught up + + // calculate percent of current segment that distance is between + double pctAlongSeg = (sampleDistance - previousDistance) / curSegmentDistance; + // add an elevation coordinate + coordList.add( + new Coordinate( + sampleDistance, + getElevation( + coverage, + new Coordinate( + x1 + (pctAlongSeg * (x2 - x1)), + y1 + (pctAlongSeg * (y2 - y1)) + ) + ) + ) + ); + sampleDistance += distanceBetweenSamplesM; + } + previousDistance = edgeLenM; + x1 = x2; + y1 = y2; } - Coordinate internal = getPointAlongEdge(coords, edgeLenM, x / edgeLenM); - coordList.add(new Coordinate(x, getElevation(internal))); - } + // remove final-segment sample if it is less than half the distance between samples + if (edgeLenM - coordList.get(coordList.size() - 1).x < distanceBetweenSamplesM / 2) { + coordList.remove(coordList.size() - 1); + } - // final sample (x = edge length) - coordList.add(new Coordinate(edgeLenM, getElevation(coords[coords.length - 1]))); + // final sample (x = edge length) + coordList.add(new Coordinate(edgeLenM, getElevation(coverage, coords[coords.length - 1]))); - // construct the PCS - Coordinate coordArr[] = new Coordinate[coordList.size()]; - PackedCoordinateSequence elevPCS = new PackedCoordinateSequence.Double( - coordList.toArray(coordArr)); + // construct the PCS + Coordinate coordArr[] = new Coordinate[coordList.size()]; + PackedCoordinateSequence elevPCS = new PackedCoordinateSequence.Double( + coordList.toArray(coordArr)); - if(ee.setElevationProfile(elevPCS, false)) { - log.trace(graph.addBuilderAnnotation(new ElevationFlattened(ee))); + setEdgeElevationProfile(ee, elevPCS, graph); + } catch (PointOutsideCoverageException | TransformException e) { + log.debug("Error processing elevation for edge: {} due to error: {}", ee, e); } } - /** - * Returns a coordinate along a path located at a specific point indicated by the percentage of - * distance covered from start to end. - * - * @param coords the list of (x,y) coordinates that form the path - * @param length the total length of the path - * @param t the percentage (ranges from 0 to 1) - * @return the (x,y) coordinate at t - */ - public Coordinate getPointAlongEdge(Coordinate[] coords, double length, double t) { - - double pctThrough = 0; // current percentage of the edge length traversed - - // endpoints of current segment within edge: - double x1 = coords[0].x, y1 = coords[0].y, x2, y2; - - for (int i = 1; i < coords.length - 1; i++) { // loop through inner points - Coordinate innerPt = coords[i]; - x2 = innerPt.x; - y2 = innerPt.y; - - // percentage of total edge length represented by current segment: - double pct = SphericalDistanceLibrary.distance(y1, x1, y2, x2) / length; - - if (pctThrough + pct > t) { // if current segment contains 't,' we're done - double pctAlongSeg = (t - pctThrough) / pct; - return new Coordinate(x1 + (pctAlongSeg * (x2 - x1)), y1 - + (pctAlongSeg * (y2 - y1))); + private void setEdgeElevationProfile(StreetWithElevationEdge ee, PackedCoordinateSequence elevPCS, Graph graph) { + if(ee.setElevationProfile(elevPCS, false)) { + synchronized (graph) { + log.trace(graph.addBuilderAnnotation(new ElevationFlattened(ee))); } - - pctThrough += pct; - x1 = x2; - y1 = y2; } - - // handle the final segment separately - x2 = coords[coords.length - 1].x; - y2 = coords[coords.length - 1].y; - - double pct = SphericalDistanceLibrary.distance(y1, x1, y2, x2) / length; - double pctAlongSeg = (t - pctThrough) / pct; - - return new Coordinate(x1 + (pctAlongSeg * (x2 - x1)), y1 + (pctAlongSeg * (y2 - y1))); } /** * Method for retrieving the elevation at a given Coordinate. - * + * + * @param coverage the specific Coverage instance to use in order to avoid competition between threads * @param c the coordinate (NAD83) * @return elevation in meters */ - private double getElevation(Coordinate c) { - return getElevation(c.x, c.y); + private double getElevation(Coverage coverage, Coordinate c) throws PointOutsideCoverageException, TransformException { + return getElevation(coverage, c.x, c.y); } /** * Method for retrieving the elevation at a given (x, y) pair. - * + * + * @param coverage the specific Coverage instance to use in order to avoid competition between threads * @param x the query longitude (NAD83) * @param y the query latitude (NAD83) * @return elevation in meters */ - private double getElevation(double x, double y) { + private double getElevation(Coverage coverage, double x, double y) throws PointOutsideCoverageException, TransformException { double values[] = new double[1]; try { // We specify a CRS here because otherwise the coordinates are assumed to be in the coverage's native CRS. @@ -464,16 +603,58 @@ private double getElevation(double x, double y) { // for DefaultGeographicCRS.WGS84, but OTP is using (long, lat) throughout and assumes unprojected DEM // rasters to also use (long, lat). coverage.evaluate(new DirectPosition2D(GeometryUtils.WGS84_XY, x, y), values); - } catch (org.opengis.coverage.PointOutsideCoverageException e) { - nPointsOutsideDEM += 1; + } catch (PointOutsideCoverageException e) { + nPointsOutsideDEM.incrementAndGet(); + throw e; } - nPointsEvaluated += 1; - return values[0]; + nPointsEvaluated.incrementAndGet(); + return values[0] - (includeEllipsoidToGeoidDifference ? getApproximateEllipsoidToGeoidDifference(y, x) : 0); + } + + /** + * The Calculation of the EllipsoidToGeoidDifference is a very expensive operation, so the resulting values are + * cached based on the coordinate values up to 2 significant digits. Two significant digits are often more than enough + * for most parts of the world, but is useful for certain areas that have dramatic changes. Since the values are + * computed once and cached, it has almost no affect on performance to have this level of detail. + * See this image for an approximate mapping of these difference values: + * https://earth-info.nga.mil/GandG/images/ww15mgh2.gif + * + * @param y latitude + * @param x longitude + */ + private double getApproximateEllipsoidToGeoidDifference(double y, double x) throws TransformException { + int geoidDifferenceCoordinateValueMultiplier = 100; + int xVal = (int) Math.round(x * geoidDifferenceCoordinateValueMultiplier); + int yVal = (int) Math.round(y * geoidDifferenceCoordinateValueMultiplier); + // create a hash value that can be used to look up the value for the given rounded coordinate. The expected + // value of xVal should never be less than -18000 (-180 * 100) or more than 18000 (180 * 100), so multiply the + // yVal by a prime number of a magnitude larger so that there won't be any hash collisions. + int hash = yVal * 104729 + xVal; + Double difference = geoidDifferenceCache.get(hash); + if (difference == null) { + difference = computeEllipsoidToGeoidDifference( + yVal / (1.0 * geoidDifferenceCoordinateValueMultiplier), + xVal / (1.0 * geoidDifferenceCoordinateValueMultiplier) + ); + geoidDifferenceCache.put(hash, difference); + } + return difference; } @Override public void checkInputs() { gridCoverageFactory.checkInputs(); + + // check for the existence of cached elevation data. + if (readCachedElevations) { + if (Files.exists(cachedElevationsFile.toPath())) { + log.info("Cached elevations file found!"); + } else { + log.warn("No cached elevations file found or read access not allowed! Unable to load in cached elevations. This could take a while..."); + } + } else { + log.warn("Not using cached elevations! This could take a while..."); + } } } \ No newline at end of file diff --git a/src/main/java/org/opentripplanner/graph_builder/module/ned/GeotiffGridCoverageFactoryImpl.java b/src/main/java/org/opentripplanner/graph_builder/module/ned/GeotiffGridCoverageFactoryImpl.java index 5284588402a..e25c4f006cd 100644 --- a/src/main/java/org/opentripplanner/graph_builder/module/ned/GeotiffGridCoverageFactoryImpl.java +++ b/src/main/java/org/opentripplanner/graph_builder/module/ned/GeotiffGridCoverageFactoryImpl.java @@ -1,6 +1,7 @@ package org.opentripplanner.graph_builder.module.ned; import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.Interpolator2D; import org.geotools.factory.Hints; import org.geotools.gce.geotiff.GeoTiffFormat; import org.geotools.gce.geotiff.GeoTiffReader; @@ -9,6 +10,7 @@ import org.slf4j.Logger; import org.slf4j.LoggerFactory; +import javax.media.jai.InterpolationBilinear; import java.io.File; import java.io.IOException; @@ -20,7 +22,6 @@ public class GeotiffGridCoverageFactoryImpl implements ElevationGridCoverageFact private static final Logger LOG = LoggerFactory.getLogger(GeotiffGridCoverageFactoryImpl.class); private final File path; - private GridCoverage2D coverage; public GeotiffGridCoverageFactoryImpl(File path) { this.path = path; @@ -28,6 +29,7 @@ public GeotiffGridCoverageFactoryImpl(File path) { @Override public GridCoverage2D getGridCoverage() { + GridCoverage2D coverage; try { // There is a serious standardization failure around the axis order of WGS84. See issue #1930. // GeoTools assumes strict EPSG axis order of (latitude, longitude) unless told otherwise. @@ -37,7 +39,9 @@ public GridCoverage2D getGridCoverage() { GeoTiffFormat format = new GeoTiffFormat(); GeoTiffReader reader = format.getReader(path, forceLongLat); coverage = reader.read(null); - LOG.info("Elevation model CRS is: {}", coverage.getCoordinateReferenceSystem2D()); + LOG.debug("Elevation model CRS is: {}", coverage.getCoordinateReferenceSystem2D()); + // TODO might bicubic interpolation give better results? + coverage = Interpolator2D.create(coverage, new InterpolationBilinear()); } catch (IOException e) { throw new RuntimeException("Error getting coverage automatically. ", e); } @@ -52,7 +56,7 @@ public void checkInputs() { } @Override - public void setGraph(Graph graph) { + public void fetchData(Graph graph) { //nothing to do here } diff --git a/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDDownloader.java b/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDDownloader.java index fccc27bb1b3..e7f4a6b066d 100644 --- a/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDDownloader.java +++ b/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDDownloader.java @@ -45,15 +45,7 @@ public class NEDDownloader implements NEDTileSource { private double _lonXStep = 0.16; - @Override - public void setGraph(Graph graph) { - this.graph = graph; - } - - @Override - public void setCacheDirectory(File cacheDirectory) { - this.cacheDirectory = cacheDirectory; - } + private List nedTiles; private File getPathToNEDArchive(String key) { if (!cacheDirectory.exists()) { @@ -294,7 +286,15 @@ private String getKey(URL url) { @Override public List getNEDTiles() { + return nedTiles; + } + + @Override public void fetchData(Graph graph, File cacheDirectory) { + this.graph = graph; + this.cacheDirectory = cacheDirectory; + log.info("Downloading NED elevation data (or fetching it from local cache)."); + List urls = getDownloadURLsCached(); List files = new ArrayList(); int tileCount = 0; @@ -336,7 +336,7 @@ public List getNEDTiles() { } log.error("Unable to download a NED tile after 5 requests."); } - return files; + nedTiles = files; } private List getDownloadURLsCached() { diff --git a/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDGridCoverageFactoryImpl.java b/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDGridCoverageFactoryImpl.java index 912c7218d11..693ec5ed6ec 100644 --- a/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDGridCoverageFactoryImpl.java +++ b/src/main/java/org/opentripplanner/graph_builder/module/ned/NEDGridCoverageFactoryImpl.java @@ -29,19 +29,20 @@ public class NEDGridCoverageFactoryImpl implements ElevationGridCoverageFactory private static final Logger LOG = LoggerFactory.getLogger(NEDGridCoverageFactoryImpl.class); - private Graph graph; + private final File cacheDirectory; - /** All tiles for the DEM stitched into a single coverage. */ - UnifiedGridCoverage unifiedCoverage = null; - - private File cacheDirectory; - - public NEDTileSource tileSource = new NEDDownloader(); + public final NEDTileSource tileSource; private List datums; public NEDGridCoverageFactoryImpl(File cacheDirectory) { this.cacheDirectory = cacheDirectory; + this.tileSource = new NEDDownloader(); + } + + public NEDGridCoverageFactoryImpl(File cacheDirectory, NEDTileSource tileSource) { + this.cacheDirectory = cacheDirectory; + this.tileSource = tileSource; } private static final String[] DATUM_FILENAMES = {"g2012a00.gtx", "g2012g00.gtx", "g2012h00.gtx", "g2012p00.gtx", "g2012s00.gtx", "g2012u00.gtx"}; @@ -101,24 +102,19 @@ private void loadVerticalDatum () { } } - /** @return a GeoTools grid coverage for the entire area of interest, lazy-creating it on the first call. */ + /** @return a GeoTools grid coverage for the entire area of interest. */ public Coverage getGridCoverage() { - if (unifiedCoverage == null) { - loadVerticalDatum(); - tileSource.setGraph(graph); - tileSource.setCacheDirectory(cacheDirectory); - List paths = tileSource.getNEDTiles(); - // Make one grid coverage for each NED tile, adding them all to a single UnifiedGridCoverage. - for (File path : paths) { - GeotiffGridCoverageFactoryImpl factory = new GeotiffGridCoverageFactoryImpl(path); - // TODO might bicubic interpolation give better results? - GridCoverage2D regionCoverage = Interpolator2D.create(factory.getGridCoverage(), - new InterpolationBilinear()); - if (unifiedCoverage == null) { - unifiedCoverage = new UnifiedGridCoverage("unified", regionCoverage, datums); - } else { - unifiedCoverage.add(regionCoverage); - } + loadVerticalDatum(); + List paths = tileSource.getNEDTiles(); + UnifiedGridCoverage unifiedCoverage = null; + // Make one grid coverage for each NED tile, adding them all to a single UnifiedGridCoverage. + for (File path : paths) { + GeotiffGridCoverageFactoryImpl factory = new GeotiffGridCoverageFactoryImpl(path); + GridCoverage2D regionCoverage = factory.getGridCoverage(); + if (unifiedCoverage == null) { + unifiedCoverage = new UnifiedGridCoverage("unified", regionCoverage, datums); + } else { + unifiedCoverage.add(regionCoverage); } } return unifiedCoverage; @@ -181,9 +177,12 @@ public void checkInputs() { } - /** Set the graph that will be used to determine the extent of the NED. */ + /** + * Verify that the needed elevation data exists in the cache and if it does not exist, try to download it. The graph + * is used to determine the extent of the NED. + */ @Override - public void setGraph(Graph graph) { - this.graph = graph; + public void fetchData(Graph graph) { + tileSource.fetchData(graph, cacheDirectory); } } \ No newline at end of file diff --git a/src/main/java/org/opentripplanner/graph_builder/module/ned/UnifiedGridCoverage.java b/src/main/java/org/opentripplanner/graph_builder/module/ned/UnifiedGridCoverage.java index b1ee0201a64..0a14b47b76a 100644 --- a/src/main/java/org/opentripplanner/graph_builder/module/ned/UnifiedGridCoverage.java +++ b/src/main/java/org/opentripplanner/graph_builder/module/ned/UnifiedGridCoverage.java @@ -1,8 +1,12 @@ package org.opentripplanner.graph_builder.module.ned; +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; +import com.vividsolutions.jts.index.SpatialIndex; +import com.vividsolutions.jts.index.strtree.STRtree; import org.geotools.coverage.AbstractCoverage; import org.geotools.coverage.grid.GridCoverage2D; -import org.geotools.geometry.GeneralEnvelope; +import org.geotools.geometry.jts.ReferencedEnvelope; import org.opengis.coverage.CannotEvaluateException; import org.opengis.coverage.Coverage; import org.opengis.coverage.PointOutsideCoverageException; @@ -26,9 +30,16 @@ public class UnifiedGridCoverage extends AbstractCoverage { private static final long serialVersionUID = -7798801307087575896L; private static Logger log = LoggerFactory.getLogger(UnifiedGridCoverage.class); - - private ArrayList regions; + /** + * A spatial index of the intersection of all regions and datums. For smaller-scale deployments, this spatial index + * might perform more slowly than manually iterating over each region and datum. However, in larger regions, this + * can result in 20+% more calculations being done during the same time compared to the brute-force method. In + * smaller-scale deployments since the overall time is not as long, we leave this in here for the benefit of larger + * regions where this will result in much better performance. + */ + private final SpatialIndex datumRegionIndex; + private ArrayList regions; private List datums; /** @@ -36,52 +47,38 @@ public class UnifiedGridCoverage extends AbstractCoverage { * in the same way. However, the superclass constructor (AbstractCoverage) needs a coverage to copy properties from. * So the first sub-coverage needs to be passed in at construction time. */ - protected UnifiedGridCoverage(CharSequence name, Coverage coverage, List datums) { + protected UnifiedGridCoverage(CharSequence name, GridCoverage2D coverage, List datums) { super(name, coverage); - regions = new ArrayList(); - regions.add(coverage); + regions = new ArrayList<>(); this.datums = datums; + datumRegionIndex = new STRtree(); + // Add first coverage to list of regions/spatial index. + this.add(coverage); } @Override - public Object evaluate(DirectPosition point) throws PointOutsideCoverageException, CannotEvaluateException { + public Object evaluate(DirectPosition point) throws CannotEvaluateException { /* we don't use this function, we use evaluate(DirectPosition point, double[] values) */ return null; } - public double[] evaluate(DirectPosition point, double[] values) - throws PointOutsideCoverageException, CannotEvaluateException { - - for (Coverage region : regions) { - // GeneralEnvelope has a contains method, OpenGIS Envelope does not - GeneralEnvelope env = ((GeneralEnvelope)region.getEnvelope()); - // Check envelope to avoid incurring exception construction overhead (PointOutsideCoverageException), - // especially important when there are many regions. - if (env.contains(point)) { - double[] result; - double x = point.getOrdinate(0); - double y = point.getOrdinate(1); - try { - result = region.evaluate(point, values); - // TODO It might be faster to put all the datums and Coverage regions into a spatial index instead of iterating. - for (VerticalDatum datum : datums) { - if (datum.covers(x, y)) { - result[0] += datum.interpolatedHeight(x, y); - return result; - } - } - //if we get here, all vdatums failed. - log.error("Failed to convert elevation at " + y + ", " + x + " from NAVD88 to NAD83"); - } catch (PointOutsideCoverageException e) { - continue; - } - return result; - } + /** + * Calculate the elevation at a given point + */ + public double[] evaluate(DirectPosition point, double[] values) throws CannotEvaluateException { + double x = point.getOrdinate(0); + double y = point.getOrdinate(1); + Coordinate pointCoordinate = new Coordinate(x, y); + Envelope envelope = new Envelope(pointCoordinate); + List coverageCandidates = datumRegionIndex.query(envelope); + if (coverageCandidates.size() > 0) { + // Found a match for coverage/datum. + DatumRegion datumRegion = coverageCandidates.get(0); + double[] result = datumRegion.region.evaluate(point, values); + result[0] += datumRegion.datum.interpolatedHeight(x, y); + return result; } - /* not found */ - log.warn("Point not found: " + point); - - return null; + throw new PointOutsideCoverageException("Point not found: " + point); } @Override @@ -95,7 +92,24 @@ public SampleDimension getSampleDimension(int index) throws IndexOutOfBoundsExce } public void add(GridCoverage2D regionCoverage) { + // Iterate over datums to find intersection envelope with each region and add to spatial index. + for (VerticalDatum datum : datums) { + Envelope datumEnvelope = new Envelope(datum.lowerLeftLongitude, datum.lowerLeftLongitude + datum.deltaLongitude, datum.lowerLeftLatitude, datum.lowerLeftLatitude + datum.deltaLatitude); + ReferencedEnvelope regionEnvelope = new ReferencedEnvelope(regionCoverage.getEnvelope()); + Envelope intersection = regionEnvelope.intersection(datumEnvelope); + datumRegionIndex.insert(intersection, new DatumRegion(datum, regionCoverage)); + } regions.add(regionCoverage); } + public class DatumRegion { + public final VerticalDatum datum; + public final Coverage region; + + public DatumRegion (VerticalDatum datum, Coverage region) { + this.datum = datum; + this.region = region; + } + } + } diff --git a/src/main/java/org/opentripplanner/graph_builder/services/ned/ElevationGridCoverageFactory.java b/src/main/java/org/opentripplanner/graph_builder/services/ned/ElevationGridCoverageFactory.java index 25ab2219b4e..e70db6aaa19 100644 --- a/src/main/java/org/opentripplanner/graph_builder/services/ned/ElevationGridCoverageFactory.java +++ b/src/main/java/org/opentripplanner/graph_builder/services/ned/ElevationGridCoverageFactory.java @@ -13,9 +13,11 @@ */ public interface ElevationGridCoverageFactory { + /** Creates a new coverage instance from files already fetched */ public Coverage getGridCoverage(); public void checkInputs(); - public void setGraph(Graph graph); + /** Sets the graph of the factory and initiates the fetching of data that is not present in the cache */ + public void fetchData(Graph graph); } \ No newline at end of file diff --git a/src/main/java/org/opentripplanner/graph_builder/services/ned/NEDTileSource.java b/src/main/java/org/opentripplanner/graph_builder/services/ned/NEDTileSource.java index 276d067cdbb..3f7c4cdaf54 100644 --- a/src/main/java/org/opentripplanner/graph_builder/services/ned/NEDTileSource.java +++ b/src/main/java/org/opentripplanner/graph_builder/services/ned/NEDTileSource.java @@ -12,15 +12,12 @@ */ public interface NEDTileSource { - public abstract void setGraph(Graph graph); - /** - * The cache directory stores NED tiles. It is crucial that this be somewhere permanent - * with plenty of disk space. Don't use /tmp -- the downloading process takes a long time + * Fetches all of the required data and stores it in the specified cache directory. It is crucial that this be + * somewhere permanent with plenty of disk space. Don't use /tmp -- the downloading process takes a long time * and you don't want to repeat it if at all possible. - * @param cacheDirectory */ - public abstract void setCacheDirectory(File cacheDirectory); + public abstract void fetchData(Graph graph, File cacheDirectory); /** * Download all the NED tiles into the cache. diff --git a/src/main/java/org/opentripplanner/routing/edgetype/StreetWithElevationEdge.java b/src/main/java/org/opentripplanner/routing/edgetype/StreetWithElevationEdge.java index 1efd6c80e89..af1de39c7b7 100644 --- a/src/main/java/org/opentripplanner/routing/edgetype/StreetWithElevationEdge.java +++ b/src/main/java/org/opentripplanner/routing/edgetype/StreetWithElevationEdge.java @@ -90,6 +90,8 @@ public PackedCoordinateSequence getElevationProfile() { return CompactElevationProfile.uncompactElevationProfile(packedElevationProfile); } + public boolean hasPackedElevationProfile () { return packedElevationProfile != null; } + @Override public boolean isElevationFlattened() { return flattened; diff --git a/src/main/java/org/opentripplanner/standalone/GraphBuilderParameters.java b/src/main/java/org/opentripplanner/standalone/GraphBuilderParameters.java index 1b97cf8ee14..64f112500e4 100644 --- a/src/main/java/org/opentripplanner/standalone/GraphBuilderParameters.java +++ b/src/main/java/org/opentripplanner/standalone/GraphBuilderParameters.java @@ -1,5 +1,6 @@ package org.opentripplanner.standalone; +import org.opentripplanner.api.common.RoutingResource; import org.opentripplanner.graph_builder.module.osm.WayPropertySetSource; import org.opentripplanner.graph_builder.services.osm.CustomNamer; import org.opentripplanner.profile.StopClusterMode; @@ -185,6 +186,29 @@ public class GraphBuilderParameters { */ public String micromobilityDropoffRestrictionsUrlOrFile; + /** + * When set to true (it is by default), the elevation module will attempt to read this file in order to reuse + * calculations of elevation data for various coordinate sequences instead of recalculating them all over again. + */ + public boolean readCachedElevations; + + /** + * When set to true (it is false by default), the elevation module will create a file of a lookup map of the + * LineStrings and the corresponding calculated elevation data for those coordinates. Subsequent graph builds can + * reuse the data in this file to avoid recalculating all the elevation data again. + */ + public boolean writeCachedElevations; + + /** + * When set to true (it is false by default), the elevation module will include the Ellipsoid to Geiod difference in + * the calculations of every point along every StreetWithElevationEdge in the graph. + * + * NOTE: if this is set to true for graph building, make sure to not set the value of + * {@link RoutingResource#geoidElevation} to true otherwise OTP will add this geoid value again to all of the + * elevation values in the street edges. + */ + public boolean includeEllipsoidToGeoidDifference; + /** * Set all parameters from the given Jackson JSON tree, applying defaults. * Supplying MissingNode.getInstance() will cause all the defaults to be applied. @@ -222,6 +246,9 @@ public GraphBuilderParameters(JsonNode config) { extraEdgesStopPlatformLink = config.path("extraEdgesStopPlatformLink").asBoolean(false); micromobilityTravelRestrictionsUrlOrFile = config.path("micromobilityTravelRestrictionsUrlOrFile").asText(); micromobilityDropoffRestrictionsUrlOrFile = config.path("micromobilityDropoffRestrictionsUrlOrFile").asText(); + readCachedElevations = config.path("readCachedElevations").asBoolean(true); + writeCachedElevations = config.path("writeCachedElevations").asBoolean(false); + includeEllipsoidToGeoidDifference = config.path("includeEllipsoidToGeoidDifference").asBoolean(false); } diff --git a/src/main/java/org/opentripplanner/util/ElevationUtils.java b/src/main/java/org/opentripplanner/util/ElevationUtils.java index aa3d46dcd5e..a64c9ab0bc0 100644 --- a/src/main/java/org/opentripplanner/util/ElevationUtils.java +++ b/src/main/java/org/opentripplanner/util/ElevationUtils.java @@ -13,8 +13,22 @@ public class ElevationUtils { + // Set up a MathTransform based on the EarthGravitationalModel + private static MathTransform mt; + static { + try { + mt = new DefaultMathTransformFactory().createParameterizedTransform( + new EarthGravitationalModel.Provider().getParameters().createValue() + ); + } catch (FactoryException e) { + e.printStackTrace(); + } + } + /** - * Computes the difference between the ellipsoid and geoid at a specified lat/lon using Geotools EarthGravitationalModel + * Computes the difference between the ellipsoid and geoid at a specified lat/lon using Geotools + * EarthGravitationalModel. For unknown reasons, this method can produce incorrect results if called at the same + * time from multiple threads, so the method has been made synchronized. * * @param lat * @param lon @@ -22,13 +36,7 @@ public class ElevationUtils { * @throws FactoryException * @throws TransformException */ - - public static double computeEllipsoidToGeoidDifference(double lat, double lon) throws FactoryException, TransformException { - // Set up a MathTransform based on the EarthGravitationalModel - EarthGravitationalModel.Provider provider = new EarthGravitationalModel.Provider(); - DefaultMathTransformFactory factory = new DefaultMathTransformFactory(); - MathTransform mt = factory.createParameterizedTransform(provider.getParameters().createValue()); - + public static synchronized double computeEllipsoidToGeoidDifference(double lat, double lon) throws TransformException { // Compute the offset DirectPosition3D dest = new DirectPosition3D(); mt.transform(new DirectPosition3D(lon, lat, 0), dest); diff --git a/src/test/java/org/opentripplanner/graph_builder/module/ElevationModuleTest.java b/src/test/java/org/opentripplanner/graph_builder/module/ElevationModuleTest.java new file mode 100644 index 00000000000..1cbf3d8e53c --- /dev/null +++ b/src/test/java/org/opentripplanner/graph_builder/module/ElevationModuleTest.java @@ -0,0 +1,107 @@ +package org.opentripplanner.graph_builder.module; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.LineString; +import org.junit.Ignore; +import org.junit.Test; +import org.opentripplanner.common.geometry.GeometryUtils; +import org.opentripplanner.common.geometry.SphericalDistanceLibrary; +import org.opentripplanner.graph_builder.module.ned.DegreeGridNEDTileSource; +import org.opentripplanner.graph_builder.module.ned.ElevationModule; +import org.opentripplanner.graph_builder.module.ned.NEDGridCoverageFactoryImpl; +import org.opentripplanner.routing.edgetype.StreetTraversalPermission; +import org.opentripplanner.routing.edgetype.StreetWithElevationEdge; +import org.opentripplanner.routing.graph.Graph; +import org.opentripplanner.routing.vertextype.OsmVertex; +import org.opentripplanner.routing.vertextype.StreetVertex; +import org.opentripplanner.standalone.S3BucketConfig; + +import java.io.File; + +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertNotNull; + +public class ElevationModuleTest { + + /** + * Tests whether elevation can be properly set using s3 bucket tiles. This test is ignored by default because it + * would require checking in a file to this repository that is about 450mb large. + * + * To setup this test, copy everything in the /var/otp/cache/ned into + * .../src/test/resources/org/opentripplanner/graph_builder/module/ned/ + */ + @Test + @Ignore + public void testSetElevationOnEdgesUsingS3BucketTiles() { + // create a graph with a StreetWithElevationEdge + Graph graph = new Graph(); + OsmVertex from = new OsmVertex(graph, "from", -122.6932051, 45.5122964, 40513757); + OsmVertex to = new OsmVertex(graph, "to", -122.6903532, 45.5115309, 1677595882); + LineString geometry = GeometryUtils.makeLineString( + -122.6932051, + 45.5122964, + -122.6931118, + 45.5122687, + -122.692677, + 45.512204, + -122.692443, + 45.512142, + -122.6923995, + 45.5121229, + -122.692359, + 45.512096, + -122.692294, + 45.512024, + -122.692219, + 45.511961, + -122.69212, + 45.511908, + -122.6920612, + 45.5118889, + -122.6920047, + 45.5118749, + -122.6919636, + 45.5118711, + -122.6918977, + 45.5118684, + -122.6918211, + 45.5118662, + -122.6917612, + 45.5118599, + -122.6916928, + 45.5118503, + -122.691241, + 45.511774, + -122.6903532, + 45.5115309 + ); + Coordinate[] coordinates = geometry.getCoordinates(); + double length = 0; + for (int i = 1; i < coordinates.length; ++i) { + length += SphericalDistanceLibrary.distance(coordinates[i - 1], coordinates[i]); + } + StreetWithElevationEdge edge = new StreetWithElevationEdge( + from, + to, + geometry, + "Southwest College St", + length, + StreetTraversalPermission.ALL, + false + ); + + // create the elevation module + File cacheDirectory = new File(ElevationModuleTest.class.getResource("ned").getFile()); + DegreeGridNEDTileSource awsTileSource = new DegreeGridNEDTileSource(); + NEDGridCoverageFactoryImpl gcf = new NEDGridCoverageFactoryImpl(cacheDirectory, awsTileSource); + ElevationModule elevationModule = new ElevationModule(gcf); + + // build to graph to execute the elevation module + elevationModule.checkInputs(); + elevationModule.buildGraph(graph, new GraphBuilderModuleSummary(elevationModule)); + + // verify that elevation data has been set on the StreetWithElevationEdge + assertNotNull(edge.getElevationProfile()); + assertEquals(133.95, edge.getElevationProfile().getY(1), 0.1); + } +}