Skip to content

Commit

Permalink
set longitude,latitude axis order when working with DEMs #1930
Browse files Browse the repository at this point in the history 
also provides a static copy of WGS84 in xy axis order
and logs warnings when too many elevation evaluations fail
  • Loading branch information
@abyrd
abyrd committed May 28, 2015
1 parent e35ceb7 commit 77c30fe
Show file tree
Hide file tree
Showing 5 changed files with 82 additions and 48 deletions.
28 changes: 19 additions & 9 deletions src/main/java/org/opentripplanner/common/geometry/GeometryUtils.java
View file
Expand Up @@ -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;
Expand All @@ -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];
Expand Down
25 changes: 22 additions & 3 deletions src/main/java/org/opentripplanner/graph_builder/module/ned/ElevationModule.java
View file
Expand Up @@ -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;
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -111,8 +116,16 @@ public void buildGraph(Graph graph, HashMap<Class<?>, 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.");
}
}
}
}
}
Expand Down Expand Up @@ -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];
}

Expand Down
37 changes: 17 additions & 20 deletions ...ain/java/org/opentripplanner/graph_builder/module/ned/GeotiffGridCoverageFactoryImpl.java
View file
Expand Up @@ -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;
}

Expand Down
20 changes: 11 additions & 9 deletions src/main/java/org/opentripplanner/graph_builder/module/ned/NEDGridCoverageFactoryImpl.java
View file
Expand Up @@ -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;

Expand Down Expand Up @@ -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<File> 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;
}

/**
Expand Down
20 changes: 13 additions & 7 deletions src/main/java/org/opentripplanner/graph_builder/module/ned/UnifiedGridCoverage.java
View file
Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -44,6 +44,11 @@ public class UnifiedGridCoverage extends AbstractCoverage {

private List<VerticalDatum> 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<VerticalDatum> datums) {
super(name, coverage);
regions = new ArrayList<Coverage>();
Expand All @@ -52,9 +57,8 @@ protected UnifiedGridCoverage(CharSequence name, Coverage coverage, List<Vertica
}

@Override
public Object evaluate(DirectPosition point) throws PointOutsideCoverageException,
CannotEvaluateException {
/* we don't care about this */
public Object evaluate(DirectPosition point) throws PointOutsideCoverageException, CannotEvaluateException {
/* we don't use this function, we use evaluate(DirectPosition point, double[] values) */
return null;
}

Expand All @@ -64,13 +68,15 @@ public double[] evaluate(DirectPosition point, double[] values)
for (Coverage region : regions) {
// GeneralEnvelope has a contains method, OpenGIS Envelope does not
GeneralEnvelope env = ((GeneralEnvelope)region.getEnvelope());
// avoid incurring exception construction overhead when there are many regions
// 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);
Expand Down

0 comments on commit 77c30fe

Please sign in to comment.