diff --git a/src/main/java/org/opentripplanner/common/geometry/GeometryUtils.java b/src/main/java/org/opentripplanner/common/geometry/GeometryUtils.java index e6b48bc06de..15fc598ae28 100644 --- a/src/main/java/org/opentripplanner/common/geometry/GeometryUtils.java +++ b/src/main/java/org/opentripplanner/common/geometry/GeometryUtils.java @@ -13,13 +13,6 @@ the License, or (at your option) any later version. package org.opentripplanner.common.geometry; -import java.util.List; - -import org.geojson.GeoJsonObject; -import org.geojson.LngLatAlt; -import org.opentripplanner.analyst.UnsupportedGeometryException; -import org.opentripplanner.common.model.P2; - import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.CoordinateSequence; import com.vividsolutions.jts.geom.CoordinateSequenceFactory; @@ -31,13 +24,30 @@ the License, or (at your option) any later version. import com.vividsolutions.jts.linearref.LengthLocationMap; import com.vividsolutions.jts.linearref.LinearLocation; import com.vividsolutions.jts.linearref.LocationIndexedLine; +import org.geojson.GeoJsonObject; +import org.geojson.LngLatAlt; +import org.geotools.referencing.CRS; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opentripplanner.analyst.UnsupportedGeometryException; +import org.opentripplanner.common.model.P2; + +import java.util.List; public class GeometryUtils { - private static CoordinateSequenceFactory csf = - new Serializable2DPackedCoordinateSequenceFactory(); + private static CoordinateSequenceFactory csf = new Serializable2DPackedCoordinateSequenceFactory(); private static GeometryFactory gf = new GeometryFactory(csf); + /** A shared copy of the WGS84 CRS with longitude-first axis order. */ + public static final CoordinateReferenceSystem WGS84_XY; + static { + try { + WGS84_XY = CRS.getAuthorityFactory(true).createCoordinateReferenceSystem("EPSG:4326"); + } catch (Exception ex) { + throw new RuntimeException("Could not create longitude-first WGS84 coordinate reference system."); + } + } + public static LineString makeLineString(double... coords) { GeometryFactory factory = getGeometryFactory(); Coordinate [] coordinates = new Coordinate[coords.length / 2]; 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 819b0b39bed..76ca950a148 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 @@ -19,6 +19,7 @@ the License, or (at your option) any later version. import org.geotools.coverage.grid.Interpolator2D; import org.geotools.geometry.DirectPosition2D; import org.opengis.coverage.Coverage; +import org.opentripplanner.common.geometry.GeometryUtils; import org.opentripplanner.common.geometry.PackedCoordinateSequence; import org.opentripplanner.common.geometry.SphericalDistanceLibrary; import org.opentripplanner.common.pqueue.BinHeap; @@ -60,6 +61,10 @@ public class ElevationModule implements GraphBuilderModule { 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 distance between samples in meters. Defaults to 10m, the approximate resolution of 1/3 * arc-second NED data. @@ -111,8 +116,16 @@ public void buildGraph(Graph graph, HashMap, Object> extra) { edgesWithElevation.add(edgeWithElevation); } nProcessed += 1; - if (nProcessed % 50000 == 0) + if (nProcessed % 50000 == 0) { log.info("set elevation on {}/{} edges", nProcessed, nTotal); + double failurePercentage = nPointsOutsideDEM / nPointsEvaluated * 100; + if (failurePercentage > 50) { + log.warn("Fetching elevation failed at {}/{} points ({}%)", + 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."); + } + } } } } @@ -431,10 +444,16 @@ private double getElevation(Coordinate c) { private double getElevation(double x, double y) { double values[] = new double[1]; try { - coverage.evaluate(new DirectPosition2D(x, y), values); + // We specify a CRS here because otherwise the coordinates are assumed to be in the coverage's native CRS. + // That assumption is fine when the coverage happens to be in longitude-first WGS84 but we want to support + // GeoTIFFs in various projections. Note that GeoTools defaults to strict EPSG axis ordering of (lat, long) + // 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) { - // skip this for now + nPointsOutsideDEM += 1; } + nPointsEvaluated += 1; return values[0]; } 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 482b8137fbb..0955e74e454 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 @@ -13,50 +13,47 @@ the License, or (at your option) any later version. package org.opentripplanner.graph_builder.module.ned; -import java.io.File; -import java.io.IOException; - import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.factory.Hints; import org.geotools.gce.geotiff.GeoTiffFormat; import org.geotools.gce.geotiff.GeoTiffReader; import org.opentripplanner.graph_builder.services.ned.ElevationGridCoverageFactory; import org.opentripplanner.routing.graph.Graph; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.io.File; +import java.io.IOException; /** * Implementation of ElevationGridCoverageFactory for Geotiff data. */ public class GeotiffGridCoverageFactoryImpl implements ElevationGridCoverageFactory { - private File path = null; - private GridCoverage2D coverage; - - public GeotiffGridCoverageFactoryImpl() { + private static final Logger LOG = LoggerFactory.getLogger(GeotiffGridCoverageFactoryImpl.class); - } + private final File path; + private GridCoverage2D coverage; public GeotiffGridCoverageFactoryImpl(File path) { this.path = path; } - public void setPath(File path) { - this.path = path; - } - @Override public GridCoverage2D getGridCoverage() { - GeoTiffFormat format = new GeoTiffFormat(); - GeoTiffReader reader = null; - try { - if (path == null) { - throw new RuntimeException("Path not set"); - } - reader = format.getReader(path); + // 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. + // Both NED and SRTM data use the longitude-first axis order, so OTP makes grid coverages + // for unprojected DEMs assuming coordinates are in (longitude, latitude) order. + Hints forceLongLat = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE); + GeoTiffFormat format = new GeoTiffFormat(); + GeoTiffReader reader = format.getReader(path, forceLongLat); coverage = reader.read(null); + LOG.info("Elevation model CRS is: {}", coverage.getCoordinateReferenceSystem2D()); } catch (IOException e) { throw new RuntimeException("Error getting coverage automatically. ", e); } - return coverage; } 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 5a07d801802..e4e810603a5 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 @@ -45,7 +45,7 @@ public class NEDGridCoverageFactoryImpl implements ElevationGridCoverageFactory private Graph graph; /** All tiles for the DEM stitched into a single coverage. */ - UnifiedGridCoverage coverage = null; + UnifiedGridCoverage unifiedCoverage = null; private File cacheDirectory; @@ -113,26 +113,28 @@ private void loadVerticalDatum () { } } } - + + /** @return a GeoTools grid coverage for the entire area of interest, lazy-creating it on the first call. */ public Coverage getGridCoverage() { - if (coverage == null) { + 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(); - factory.setPath(path); + GeotiffGridCoverageFactoryImpl factory = new GeotiffGridCoverageFactoryImpl(path); + // TODO might bicubic interpolation give better results? GridCoverage2D regionCoverage = Interpolator2D.create(factory.getGridCoverage(), new InterpolationBilinear()); - if (coverage == null) { - coverage = new UnifiedGridCoverage("unified", regionCoverage, datums); + if (unifiedCoverage == null) { + unifiedCoverage = new UnifiedGridCoverage("unified", regionCoverage, datums); } else { - coverage.add(regionCoverage); + unifiedCoverage.add(regionCoverage); } } } - return coverage; + return unifiedCoverage; } /** 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 3b106781cf4..e7a67cdb7c9 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 @@ -13,9 +13,6 @@ the License, or (at your option) any later version. package org.opentripplanner.graph_builder.module.ned; -import java.util.ArrayList; -import java.util.List; - import org.geotools.coverage.AbstractCoverage; import org.geotools.coverage.grid.GridCoverage2D; import org.geotools.geometry.GeneralEnvelope; @@ -27,6 +24,9 @@ the License, or (at your option) any later version. import org.slf4j.Logger; import org.slf4j.LoggerFactory; +import java.util.ArrayList; +import java.util.List; + /** * Stitches together multiple elevation maps into a single elevation map, * hackily. This is horrible, but the geotools way of doing things is @@ -44,6 +44,11 @@ public class UnifiedGridCoverage extends AbstractCoverage { private List datums; + /** + * It would be nice if we could construct this unified coverage with zero sub-coverages and add all sub-coverages + * 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) { super(name, coverage); regions = new ArrayList(); @@ -52,9 +57,8 @@ protected UnifiedGridCoverage(CharSequence name, Coverage coverage, List