Repository: graphhopper/graphhopper
Pull request #1953 (2020-03-13 10:10:57)

Improve resolution of elevation profiles with 3D Douglas Peucker and long edge sampling
Improve the resolution of elevation profiles by:

1. sampling points along long edges so we don't miss elevation changes over that segment
2. making Douglas Peucker polyline simplification elevation-aware so it keeps points that only impact the elevation profile and not the 2D shape

These options have the potential to increase geometry size, so I added 2 configuration options to control them:

1. `graph.elevation.long_edge_sampling_distance` which defines the distance in meters that it will sample points along long edges, for example you might set this to somewhere around 30m if your elevation provider has 30m resolution. It defaults to 0 which disables long edge sampling.
2. `graph.elevation.way_point_max_distance` which defines the elevation tolerance during polyline simplification, for example with a 1m lat/lon and elevation tolerance removes points that result in a <1m shift of the 3D polyline. It defaults to Double.MAX_VALUE which disables elevation-influenced simplification.

**Edge Geometry Size Increase**

I ran this code on `north-america_us_pennsylvania.pbf` with CGIAR and bilinear interpolation with different values of those parameters and saw these edge geometry size increases over what's on master now:

sampling distance ➡️<br> elevation max distance ⬇️  | 15m | 30m | 45m | 60m | 100m | 500m | off (∞)
-- | -- | -- | -- | -- | -- | -- | --
1m | 19% | 17% | 17% | 15% | 11% | 6% | 6%
1.5m | 11% | 11% | 9% | 9% | 6% | 2% | 2%
2m | 6% | 6% | 6% | 6% | 4% | 2% | 2%
3m | 2% | 2% | 2% | 2% | 2% | 0% | 0%
4m | 2% | 2% | 2% | 2% | 0% | 0% | 0%
5m | 0% | 0% | 0% | 0% | 0% | 0% | 0%
10m | 0% | 0% | 0% | 0% | 0% | 0% | 0%
50m | 0% | 0% | 0% | 0% | 0% | 0% | 0%
off (∞) | 0% | 0% | 0% | 0% | 0% | 0% | 0%



I ran the same test with mapzen/skadi provider from #1929. Results were similar, but the increase was slightly larger (26% instead of 19% for 15m/1m) probably because of the higher resolution of the underlying data.

**Import time performance**

I observed the "pass2" time during import. There was no change for different values of the elevation factor, but shorter long edge sampling distances did increase that time by 25% for 15m, 9% for 30m, 8% for 45m, and <2% for longer distances. Last time I imported the world pass2 took 110 minutes, so even a 25% increase is only about a 28 minute increase which seems negligible to me.

**Examples**

This [2km straight segment](https://graphhopper.com/maps/?point=41.7373489604297%2C-79.785794871728&point=41.7511768657977%2C-79.8039368496225&locale=en&vehicle=bike&weighting=fastest&elevation=true&turn_costs=false&use_miles=false&layer=TF%20Cycle) loses all elevation details on master
![image](https://user-images.githubusercontent.com/1480504/76612184-22630780-64f2-11ea-92aa-09f9aff49b88.png)

but with 30m edge sampling and 1m elevation resolution we get that elevation detail back:
![image](https://user-images.githubusercontent.com/1480504/76612269-4d4d5b80-64f2-11ea-8195-f6cad7c720d9.png)

Let me know if there are other examples or different values of the 2 parameters you want to see.

See further discussion [here](https://discuss.graphhopper.com/t/improving-elevation-profiles-for-short-routes/5244)

boldtrn: I gave this PR a try. I compared the results with an instance that doesn't use the resolution improvement, I used the following settings:
```
  graph.elevation.provider: skadi
  graph.elevation.interpolate: bilinear
  graph.elevation.long_edge_sampling_distance: 30
  graph.elevation.way_point_max_distance: 1
## -------------------------------------------------------------------
  graph.elevation.provider: skadi
  graph.elevation.interpolate: bilinear
  # Disable both, I debugged the creation in GraphHopper#init, Double.MAX_VALUE gets picked up correctly
  graph.elevation.long_edge_sampling_distance: Double.MAX_VALUE
  graph.elevation.way_point_max_distance: Double.MAX_VALUE

```

A major issue I saw is that while the elevation graph contains more details, the overall ascend and descend increases dramatically as well. 

For example in Baden-Würtemberg (Germany)

http://localhost:8989/maps/?point=48.037463%2C8.622894&point=49.135901%2C9.30542&locale=de-DE&vehicle=car&weighting=fastest&elevation=true&use_miles=false&layer=Omniscale

↗1514m ↘2013m => ↗2101m ↘2601m (150km)

![Screenshot from 2020-05-05 16-24-36](https://user-images.githubusercontent.com/1553525/81077269-0a55a600-8eed-11ea-9b5e-47fbbf427474.png)
![Screenshot from 2020-05-05 16-24-24](https://user-images.githubusercontent.com/1553525/81077272-0c1f6980-8eed-11ea-8f03-9f50ea3d62c7.png)

http://localhost:8989/maps/?point=48.019324%2C7.756348&point=48.341646%2C9.898682&locale=de-DE&vehicle=car&weighting=fastest&elevation=true&use_miles=false&layer=Omniscale

↗3135m ↘2834m => ↗4828m ↘4523m (300km)

http://localhost:8989/maps/?point=47.702368%2C9.541626&point=49.088258%2C8.640747&locale=de-DE&vehicle=car&weighting=fastest&elevation=true&use_miles=false&layer=Omniscale

↗3351m ↘3611m => ↗4322m ↘4583m (251km)

The reason for the difference is most probably, that I compared car routes that contain a lot of straight sections, before we might have removed many points along the polyline. Now we add extra points, which will also increase the elevation overall.

For a rather curvy mountainous segment the difference is a lot smaller, but the elevation increased:

http://localhost:8989/maps/?point=47.97803%2C7.847929&point=47.836263%2C7.937236&locale=de-DE&vehicle=car&weighting=fastest&elevation=true&use_miles=false&layer=Omniscale

↗1104m ↘682m => ↗1150m ↘727m (27km)

When you compare this segment with:
- Google Maps, shows ↗886m ↘468m
- Here WeGo, shows ↗928m ↘514m
- Komoot, shows ↗920m ↘510m
- GraphHopper Maps using CGIAR data ↗1063m ↘637m

BTW: I am not saying that the other routers are correct, but in my experience GraphHopper tends to overestimate the elevation compared to other routers, so I am a bit concerned when the elevation is increased even more. 

Not sure if you did some more research on this?
boldtrn: I tried to import Australia (same config as above, latest downloaded just now from Geofabrik) with this PR and received an Exception:

```
2020-05-06 09:01:30.647 [main] ERROR com.graphhopper.reader.osm.OSMReader - Couldn't properly add edge with osm ids:[-9223372036854775231, 253003953, 253003952, 253003951, 253003950, 253003949, 253003948, 253003947, 253003921]
java.lang.IllegalStateException: maximum not found in [0,9]
	at com.graphhopper.util.DouglasPeucker.subSimplify(DouglasPeucker.java:151)
	at com.graphhopper.util.DouglasPeucker.simplify(DouglasPeucker.java:104)
	at com.graphhopper.util.DouglasPeucker.simplify(DouglasPeucker.java:81)
	at com.graphhopper.util.DouglasPeucker.simplify(DouglasPeucker.java:77)
	at com.graphhopper.reader.osm.OSMReader.addEdge(OSMReader.java:703)
	at com.graphhopper.reader.osm.OSMReader.addOSMWay(OSMReader.java:648)
	at com.graphhopper.reader.osm.OSMReader.processWay(OSMReader.java:404)
	at com.graphhopper.reader.osm.OSMReader.writeOsmToGraph(OSMReader.java:274)
	at com.graphhopper.reader.osm.OSMReader.readGraph(OSMReader.java:150)
```

This seems to related to this [way](https://www.openstreetmap.org/way/23366648). Not sure what happened with the first OSMID `-9223372036854775231`. I didn't dig deeper so far, but will try to find out more specifics.
msbarry: @boldtrn thanks for running those tests!

> A major issue I saw is that while the elevation graph contains more details, the overall ascend and descend increases dramatically as well.

That is exactly the problem I'm trying to solve in my app with this change. My app currently samples up to 512 points along each route and I need to show this warning to users for overall ascent/descent over long routes since users have reported ascent/descent reported by GPS watches is substantially larger. 

![image](https://user-images.githubusercontent.com/1480504/81162314-f9dc1480-8f5a-11ea-8cd4-da3273e3e267.png)

The expectation with this change is that graphhopper will report larger elevation gains overall since it takes into accounts all of the bumps along the way. That will also amplify any noise/errors but from what I've seen though the skadi data has been fairly smooth.

I tried running a similar test on a shorter (4km) route to analyze the differences between providers more closely http://localhost:9898/maps/?point=48.221124%2C8.247481&point=48.237374%2C8.264079&locale=de-DE&vehicle=car&weighting=fastest&elevation=true&turn_costs=false&use_miles=false&layer=Omniscale

- local graphhopper with this change ↗435m ↘51m
- graphhopper.com without this change ↗430m ↘63m, some sudden ascents/descents in the profile, probably from crossing into a neighboring pixel without bilinear interpolation
- google maps ↗376m ↘0m *but* I have an old version of my site running that samples the same route using google maps elevation service and it reports ↗421m ↘47m so I suspect that the elevation profiels they report for routes are simplified quite a bit
- komoot ↗370m ↘10m looks similar overall, but digging into some of the local details that skadi/google elevation show and komoot misses I see a section where google, skadi, and topographic lines on komoot show elevation dropping from ~780 to ~760 (https://onthegomap.com/s/5bi93tvc) but komoot shows elevation steady at 760 +/-2m.

> BTW: I am not saying that the other routers are correct, but in my experience GraphHopper tends to overestimate the elevation compared to other routers, so I am a bit concerned when the elevation is increased even more.
> 
> Not sure if you did some more research on this?

I tested this out on bike rides and runs that I've done in northeastern US. The elevation profiles matched my expectation better, but I didn't compare exact ascent/descent amounts.  I suspect all of the major route providers are underestimating elevation change, possibly in a tradeoff to show less noise at the expense of losing real details? One source of graphhopper overestimating in the past could be that without bilinear interpolation some sampling noise was mixing in?

Another thought I had to reduce noise would be additional interpolation modes that take into account more of the local neighborhood than bilinear interpolation does to get an extra level of smoothness, but that could be decoupled from this change.
msbarry: > I tried to import Australia (same config as above, latest downloaded just now from Geofabrik) with this PR and received an Exception:
> 
> ```
> 2020-05-06 09:01:30.647 [main] ERROR com.graphhopper.reader.osm.OSMReader - Couldn't properly add edge with osm ids:[-9223372036854775231, 253003953, 253003952, 253003951, 253003950, 253003949, 253003948, 253003947, 253003921]
> java.lang.IllegalStateException: maximum not found in [0,9]
> 	at com.graphhopper.util.DouglasPeucker.subSimplify(DouglasPeucker.java:151)
> 	at com.graphhopper.util.DouglasPeucker.simplify(DouglasPeucker.java:104)
> 	at com.graphhopper.util.DouglasPeucker.simplify(DouglasPeucker.java:81)
> 	at com.graphhopper.util.DouglasPeucker.simplify(DouglasPeucker.java:77)
> 	at com.graphhopper.reader.osm.OSMReader.addEdge(OSMReader.java:703)
> 	at com.graphhopper.reader.osm.OSMReader.addOSMWay(OSMReader.java:648)
> 	at com.graphhopper.reader.osm.OSMReader.processWay(OSMReader.java:404)
> 	at com.graphhopper.reader.osm.OSMReader.writeOsmToGraph(OSMReader.java:274)
> 	at com.graphhopper.reader.osm.OSMReader.readGraph(OSMReader.java:150)
> ```
> 
> This seems to related to this [way](https://www.openstreetmap.org/way/23366648). Not sure what happened with the first OSMID `-9223372036854775231`. I didn't dig deeper so far, but will try to find out more specifics.

Good catch - looks like 3d douglas peucker doesn't handle start/end points being the same. I'll add some tests for this case and fix.
boldtrn: Awesome thanks, I can confirm, Australia is working fine now.
boldtrn: Thanks for the detailed feedback.

> The expectation with this change is that graphhopper will report larger elevation gains overall since it takes into accounts all of the bumps along the way. That will also amplify any noise/errors but from what I've seen though the skadi data has been fairly smooth.

One thing I have seen some times in the past is, if a way leads along a ridge / cliff, the elevation would be picked sometimes from on top of the ridge and sometimes from the lower area, which adds a lot of incorrect noise. I am not sure how much this changed with the bilinear smoothing.

A good example where you can see the issue quite well is the Danube valley. The road more or less follows the river, so it's pretty much flat. There is no serious elevation along this road. 

http://localhost:8989/maps/?point=48.079369%2C9.191995&point=48.065634%2C9.001451&locale=de-DE&elevation=true&vehicle=car&weighting=fastest&turn_costs=false&use_miles=false&layer=TF%20Outdoors

This PR reports ↗589m ↘555m, Komoot reports ↗170m ↘150m. Komoot is definitely closer to the truth at this point. We do have a video of the Danube valley on the Kurviger Youtube channel, there you can see that it's pretty flat (the last sequence is up a hill that is not part of the above route). In the video it's also easy to see that big parts of the route lead along very steep cliffs (e.g. at 3:00).

https://www.youtube.com/watch?v=IqjvdjtQclc

That said, with the bilinear smoothing this issue has been massively improved. This was a lot worse in the past.

> Have reported ascent/descent reported by GPS watches is substantially larger.

I have no GPS watch, so I don't have any experience here. But could it be that the GPS watches are summing up errors as well? For example if we assume that the GPS has an average error of 8m, this could sum up quite fast. If they also use barometric pressure this might introduce some issues as well.

I remember a blog post from Strava, but cannot find it right now I could only find some info that they will smooth the input data to reduce noise. But I found this post which looks interesting as well, showing how errors used to sum up from GPS sensors: http://regex.info/blog/2015-05-09/2568

In any case, I think we need two things:
- Filtering noise from the elevation calculation
- Improving the resolution

I think both are important and somewhat contrary, I think we need a bit of both :).
boldtrn: Maybe we should move the issue on how the elevation is summed up into a different issue?

By default, we could set the sampling to (just an example):
``` 
-Ddw.graphhopper.graph.elevation.long_edge_sampling_distance=100
-Ddw.graphhopper.graph.elevation.way_point_max_distance=5
``` 

Using that setting, I cannot see big differences in the overall ascend/descend, but for big gaps (>100m) in ways we would fill these gaps to not miss any serious elevation changes. Everyone who needs a higher resolution, could set it up that way.

I think we should run some more tests to find a good default value.

WDYT does that make sense?
msbarry: > Maybe we should move the issue on how the elevation is summed up into a different issue?
> 
> By default, we could set the sampling to (just an example):
> 
> ```
> -Ddw.graphhopper.graph.elevation.long_edge_sampling_distance=100
> -Ddw.graphhopper.graph.elevation.way_point_max_distance=5
> ```
> 
> Using that setting, I cannot see big differences in the overall ascend/descend, but for big gaps (>100m) in ways we would fill these gaps to not miss any serious elevation changes. Everyone who needs a higher resolution, could set it up that way.
> 
> I think we should run some more tests to find a good default value.
> 
> WDYT does that make sense?

I think that makes sense, `graph.elevation.way_point_max_distance` is the knob we have tune down sensitivity and filter-out noise, given the points you raise 5m seems like a reasonable default value, and also keeps edge geometry size increase negligible. For the Danube case, it looks like FWIW, Google maps elevation profile simplifies it down to just 7 elevation points over the 20km route.

I could see how pure bilinear interpolation might make it think that points near a slope like in the Danube case are actually on the slope, I'd be curious to experiment with other interpolation methods like bicubic that might be able to infer more accurately where a slope starts between elevation samples

![image](https://user-images.githubusercontent.com/1480504/81242841-2b48f480-8fdc-11ea-8167-8a25b61140dd.png)

@boldtrn do you want me to set the default values to sample distance=100/elevation max dist=5 in this review and open up a separate issue for smoothing elevation noise?
msbarry: @boldtrn so to sum up here it sounds like the path forward is to:

- try different edge sampling settings to find one that gives overall elevation comparable to other providers and enable that by default
- confirm that great circle interpolation on long points is doing what we want on those australia edges
- open a new issue for smoothing noise in elevation profiles/ascent/descent calculations

?

Also, the elevation profile seems to have disappeared from my local graphhopper web UI - is there something special you need to always have in config.yml to make that show up?
karussell: > Also, the elevation profile seems to have disappeared from my local graphhopper web UI - is there something special you need to always have in config.yml to make that show up?

Unfortunately this is a bug, likely related to #1952, but will have a closer look as it should be involved with the profiles changes.
boldtrn: > FWIW, Google maps simplified the elevation for the Danube profile to just 7 elevation points over the 20km route.

Are these the dots shown in the image below? How do you find this out?

> I'd be curious to experiment with other interpolation methods like bicubic that might be able to infer more accurately where a slope starts between elevation samples

Yes indeed, that could help in this situation.

> so to sum up here it sounds like the path forward is to

Yes that sounds good to me.

> enable that by default

I am not 100% sure about enabling it by default or just setting these values as default, when it's enabled, like through config-example. If we enable it by default, we should probably also enable bilinear smoothing by default? WDYT? @karussell?
msbarry: > Are these the dots shown in the image below? How do you find this out?

Those are the elevation sample locations, I used https://turfjs.org/docs/#pointGrid to overlay where they would be based on 3600 samples per degree.
msbarry: OK a couple of updates here:
- I got the UI missing elevation issue fixed by building assets locally and linking them in, thanks @karussell !
- great circle interpolation appears to work on those Australia segments, although paths do get off from the satellite view of the roads in the middle and likely always will unless someone drops a few intermediate points
- I tried bicubic interpolation, it reports larger total ascent/descent over bilinear, probably because the splines overshoot elevation samples - so I don't think that's a solution to help out here.
- I tried a bunch of variations on sample distance and elevation max distance. For a default value I'd probably stick to <100m sampling so we are getting more data and rely more on the elevation tolerance (3-10m) to smooth out noise.

![image](https://user-images.githubusercontent.com/1480504/81401978-79eeaf80-90fe-11ea-8eba-f0ba8e47d7ec.png)
(green/red are difference relative to baseline CGIAR without interpolation)

Also, digging into some of the provider discrepancies issues I'm seeing cases where every contour map shows the road going up a hill, but it appears the road is actually flat there (and komoot/google don't show elevation change). 

<img width="1193" alt="Screen Shot 2020-05-07 at 4 28 20 PM" src="https://user-images.githubusercontent.com/1480504/81398742-d00c2480-90f7-11ea-8637-b59806c9b7b8.png">

I think the only ways to correct for this would either be to post-process the way geometries to smooth out elevation or if you had enough user-provided GPS data snapped to OSM ways, you could use that to fix issues in elevation like in the article @boldtrn shared.  Only the former would be practical for the open source project. It looks like the `graph.elevation.smoothing` option @boldtrn added in #1220 reduces elevation changes that are reported, but it introduces some sharp ascent/descent artifacts to the profile that would be not work at least for my use-case. Do you think the new issue should be around improving that smoothing code? or a completely different approach?
otbutz: Are there any side effects for the existing elevation interpolation used for bridges and tunnels?
otbutz: Maybe we should avoid smoothing the elevation for all points in GraphElevationSmoothing and only use this approach if a certain slope threshold is exceeded?
karussell: Good catch with the GraphElevationSmoothing (@boldtrn will know more) but EdgeElevationInterpolator (bridges/tunnels) is unrelated here. For this we only want to "get rid" of the elevation data.
msbarry: I submitted #2028 to track improved elevation profile smoothing.  Let me know what else you'd like me to change on this review before merging.
boldtrn: > great circle interpolation appears to work on those Australia segments, although paths do get off from the satellite view of the roads in the middle and likely always will unless someone drops a few intermediate points. Let me know what you think about the implementation.

Yes, I think this is the expected behavior. So no worries here :).

> I tried bicubic interpolation, it reports larger total ascent/descent over bilinear, probably because the splines overshoot elevation samples - so I don't think that's a solution to help out here.

That's a pity, but thanks for trying. I was hoping this might help in this situation.

> Also, digging into some of the provider discrepancies issues I'm seeing cases where every contour map shows the road going up a hill, but it appears the road is actually flat there (and komoot/google don't show elevation change).

That is actually quite an interesting find. Do you know why this could happen? Maybe the contour lines are approximated already? Therefore, it might not be possible to fully trust them?

Thanks for creating the new issue :).
otbutz: > That is actually quite an interesting find. Do you know why this could happen? Maybe the contour lines are approximated already? Therefore, it might not be possible to fully trust them?

Might be caused by the interpolation in the elevation provider. @msbarry you could try to disable it for your tests
msbarry: > How about we go with:
> 
> ```
>   graph.elevation.long_edge_sampling_distance: 60
>   graph.elevation.way_point_max_distance: 10
> ```
> 
> These would be configurable anyway, but results for these look good straight away if anyone wants to activate this feature. WDYT?

Sounds good to me, I set those default values in the comments in `config-example.yml`.

> > That is actually quite an interesting find. Do you know why this could happen? Maybe the contour lines are approximated already? Therefore, it might not be possible to fully trust them?
> 
> Might be caused by the interpolation in the elevation provider. @msbarry you could try to disable it for your tests

I don't think it's an issue sampling the actual DEM file since this based on where the samples are located, there should be enough resolution to separate the river from the hill around that curve
![image](https://user-images.githubusercontent.com/1480504/82052518-c0598680-9689-11ea-8b2d-84cf399706fb.png)

It's probably an irregularity with the DEM data provider. According to https://github.com/tilezen/joerd/blob/master/docs/data-sources.md, skadi uses [EU-DEM](https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem) for most of Europe which is "a hybrid product based on SRTM and ASTER GDEM data fused by a weighted averaging approach." So it could be artifacts from either in the SRTM/ASTER -> EU-DEM conversion or EU-DEM to skadi conversion.
msbarry: Thanks for the review @boldtrn! @karussell any more changes you want me to make here?
msbarry: Congrats on the 1.0 release! Think this might be safe to merge soon @karussell ?
karussell: Thanks a lot @msbarry for your work here! Also thanks @boldtrn for the review!

As I understand it correctly an existing application is not affected or will it be affected if it already uses interpolation?
karussell: @otbutz 

> Are there any side effects for the existing elevation interpolation used for bridges and tunnels?

I don't think so. Basically what we do for bridges and tunnels is to get rid of elevation.

> Edit: if i'm not mistaken, we interpolate elevation at three different locations now?

Yes, but at least the bridge/tunnel part is a different purpose.
msbarry: > As I understand it correctly an existing application is not affected or will it be affected if it already uses interpolation?

Yes, existing applications should not be affected whether they use bilinear interpolation or not since the 2 new config parameters default to Double.MAX_VALUE which retains existing behavior.
msbarry: > Now I also better understand your images in [#1953 (comment)](https://github.com/graphhopper/graphhopper/pull/1953#issuecomment-625774030), i.e. red is the only area where the incline and decline increase ("got worse"). So this PR with the proposed settings in config-example.yml is always an improvement and should decrease incline&decline? (over CGIAR and interpolation only)

Those charts show increase/decrease relative to CGIAR without interpolation. 10m/60m plus bilinear interpolation proposed in config-example.yml should report an "improvement" with less total elevation gain/loss than people are used to with CGIAR and no interpolation today.

Patch:
diff --git a/api/src/main/java/com/graphhopper/util/DistanceCalc.java b/api/src/main/java/com/graphhopper/util/DistanceCalc.java
index 03f753bd0..2b9d3e208 100644
--- a/api/src/main/java/com/graphhopper/util/DistanceCalc.java
+++ b/api/src/main/java/com/graphhopper/util/DistanceCalc.java
@@ -113,6 +113,12 @@ public interface DistanceCalc {
     GHPoint projectCoordinate(double lat, double lon,
                               double distanceInMeter, double headingClockwiseFromNorth);
 
+    /**
+     * This methods creates a point (lat, lon in degrees) a fraction of the distance along the path from (lat1, lon1)
+     * to (lat2, lon2).
+     */
+    GHPoint intermediatePoint(double f, double lat1, double lon1, double lat2, double lon2);
+
     /*
      * Simple heuristic to detect if the specified two points are crossing the boundary +-180°. See
      * #667
diff --git a/api/src/main/java/com/graphhopper/util/DistanceCalcEarth.java b/api/src/main/java/com/graphhopper/util/DistanceCalcEarth.java
index 5262a5dc7..2df8be320 100644
--- a/api/src/main/java/com/graphhopper/util/DistanceCalcEarth.java
+++ b/api/src/main/java/com/graphhopper/util/DistanceCalcEarth.java
@@ -119,91 +119,71 @@ public class DistanceCalcEarth implements DistanceCalc {
     public double calcNormalizedEdgeDistance(double r_lat_deg, double r_lon_deg,
                                              double a_lat_deg, double a_lon_deg,
                                              double b_lat_deg, double b_lon_deg) {
-        return calcNormalizedEdgeDistanceNew(r_lat_deg, r_lon_deg, a_lat_deg, a_lon_deg, b_lat_deg, b_lon_deg, false);
-    }
-
-    @Override
-    public double calcNormalizedEdgeDistance3D(double r_lat_deg, double r_lon_deg, double r_ele_m,
-                                               double a_lat_deg, double a_lon_deg, double a_ele_m,
-                                               double b_lat_deg, double b_lon_deg, double b_ele_m) {
-        if (Double.isNaN(r_ele_m) || Double.isNaN(a_ele_m) || Double.isNaN(b_ele_m))
-            return calcNormalizedEdgeDistance(r_lat_deg, r_lon_deg, a_lat_deg, a_lon_deg, b_lat_deg, b_lon_deg);
-
         double shrinkFactor = calcShrinkFactor(a_lat_deg, b_lat_deg);
 
         double a_lat = a_lat_deg;
         double a_lon = a_lon_deg * shrinkFactor;
-        double a_ele = a_ele_m / METERS_PER_DEGREE;
 
         double b_lat = b_lat_deg;
         double b_lon = b_lon_deg * shrinkFactor;
-        double b_ele = b_ele_m / METERS_PER_DEGREE;
 
         double r_lat = r_lat_deg;
         double r_lon = r_lon_deg * shrinkFactor;
-        double r_ele = r_ele_m / METERS_PER_DEGREE;
 
         double delta_lon = b_lon - a_lon;
         double delta_lat = b_lat - a_lat;
-        double delta_ele = b_ele - a_ele;
 
-        double norm = delta_lon * delta_lon + delta_lat * delta_lat + delta_ele * delta_ele;
-        double factor = ((r_lon - a_lon) * delta_lon + (r_lat - a_lat) * delta_lat + (r_ele - a_ele) * delta_ele) / norm;
+        if (delta_lat == 0)
+            // special case: horizontal edge
+            return calcNormalizedDist(a_lat_deg, r_lon_deg, r_lat_deg, r_lon_deg);
 
-        // x,y,z is projection of r onto segment a-b
+        if (delta_lon == 0)
+            // special case: vertical edge
+            return calcNormalizedDist(r_lat_deg, a_lon_deg, r_lat_deg, r_lon_deg);
+
+        double norm = delta_lon * delta_lon + delta_lat * delta_lat;
+        double factor = ((r_lon - a_lon) * delta_lon + (r_lat - a_lat) * delta_lat) / norm;
+
+        // x,y is projection of r onto segment a-b
         double c_lon = a_lon + factor * delta_lon;
         double c_lat = a_lat + factor * delta_lat;
-        double c_ele_m = (a_ele + factor * delta_ele) * METERS_PER_DEGREE;
-        return calcNormalizedDist(c_lat, c_lon / shrinkFactor, r_lat_deg, r_lon_deg) + calcNormalizedDist(r_ele_m - c_ele_m);
+        return calcNormalizedDist(c_lat, c_lon / shrinkFactor, r_lat_deg, r_lon_deg);
     }
 
-    /**
-     * New edge distance calculation where no validEdgeDistance check would be necessary
-     * <p>
-     *
-     * @return the normalized distance of the query point "r" to the project point "c" onto the line
-     * segment a-b
-     */
-    public double calcNormalizedEdgeDistanceNew(double r_lat_deg, double r_lon_deg,
-                                                double a_lat_deg, double a_lon_deg,
-                                                double b_lat_deg, double b_lon_deg,
-                                                boolean reduceToSegment) {
+    @Override
+    public double calcNormalizedEdgeDistance3D(double r_lat_deg, double r_lon_deg, double r_ele_m,
+                                               double a_lat_deg, double a_lon_deg, double a_ele_m,
+                                               double b_lat_deg, double b_lon_deg, double b_ele_m) {
+        if (Double.isNaN(r_ele_m) || Double.isNaN(a_ele_m) || Double.isNaN(b_ele_m))
+            return calcNormalizedEdgeDistance(r_lat_deg, r_lon_deg, a_lat_deg, a_lon_deg, b_lat_deg, b_lon_deg);
+
         double shrinkFactor = calcShrinkFactor(a_lat_deg, b_lat_deg);
 
         double a_lat = a_lat_deg;
         double a_lon = a_lon_deg * shrinkFactor;
+        double a_ele = a_ele_m / METERS_PER_DEGREE;
 
         double b_lat = b_lat_deg;
         double b_lon = b_lon_deg * shrinkFactor;
+        double b_ele = b_ele_m / METERS_PER_DEGREE;
 
         double r_lat = r_lat_deg;
         double r_lon = r_lon_deg * shrinkFactor;
+        double r_ele = r_ele_m / METERS_PER_DEGREE;
 
         double delta_lon = b_lon - a_lon;
         double delta_lat = b_lat - a_lat;
+        double delta_ele = b_ele - a_ele;
 
-        if (delta_lat == 0)
-            // special case: horizontal edge
-            return calcNormalizedDist(a_lat_deg, r_lon_deg, r_lat_deg, r_lon_deg);
-
-        if (delta_lon == 0)
-            // special case: vertical edge
-            return calcNormalizedDist(r_lat_deg, a_lon_deg, r_lat_deg, r_lon_deg);
-
-        double norm = delta_lon * delta_lon + delta_lat * delta_lat;
-        double factor = ((r_lon - a_lon) * delta_lon + (r_lat - a_lat) * delta_lat) / norm;
+        double norm = delta_lon * delta_lon + delta_lat * delta_lat + delta_ele * delta_ele;
+        double factor = ((r_lon - a_lon) * delta_lon + (r_lat - a_lat) * delta_lat + (r_ele - a_ele) * delta_ele) / norm;
+        if (Double.isNaN(factor)) factor = 0;
 
-        // make new calculation compatible to old
-        if (reduceToSegment) {
-            if (factor > 1)
-                factor = 1;
-            else if (factor < 0)
-                factor = 0;
-        }
-        // x,y is projection of r onto segment a-b
+        // x,y,z is projection of r onto segment a-b
         double c_lon = a_lon + factor * delta_lon;
         double c_lat = a_lat + factor * delta_lat;
-        return calcNormalizedDist(c_lat, c_lon / shrinkFactor, r_lat_deg, r_lon_deg);
+        double c_ele_m = (a_ele + factor * delta_ele) * METERS_PER_DEGREE;
+        return calcNormalizedDist(c_lat, c_lon / shrinkFactor, r_lat_deg, r_lon_deg) + calcNormalizedDist(r_ele_m - c_ele_m);
     }
 
     double calcShrinkFactor(double a_lat_deg, double b_lat_deg) {
@@ -300,6 +280,41 @@ public class DistanceCalcEarth implements DistanceCalc {
         return new GHPoint(projectedLat, projectedLon);
     }
 
+    @Override
+    public GHPoint intermediatePoint(double f, double lat1, double lon1, double lat2, double lon2) {
+        double lat1radians = Math.toRadians(lat1);
+        double lon1radians = Math.toRadians(lon1);
+        double lat2radians = Math.toRadians(lat2);
+        double lon2radians = Math.toRadians(lon2);
+
+        // This formula is taken from: (http://www.movable-type.co.uk/scripts/latlong.html -> https://github.com/chrisveness/geodesy MIT)
+
+        double deltaLat = lat2radians - lat1radians;
+        double deltaLon = lon2radians - lon1radians;
+        double cosLat1 = cos(lat1radians);
+        double cosLat2 = cos(lat2radians);
+        double sinHalfDeltaLat = sin(deltaLat / 2);
+        double sinHalfDeltaLon = sin(deltaLon / 2);
+
+        double a = sinHalfDeltaLat * sinHalfDeltaLat + cosLat1 * cosLat2 * sinHalfDeltaLon * sinHalfDeltaLon;
+        double angularDistance =  2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
+        double sinDistance = sin(angularDistance);
+
+        if (angularDistance == 0) return new GHPoint(lat1, lon1);
+
+        double A = Math.sin((1-f)*angularDistance) / sinDistance;
+        double B = Math.sin(f*angularDistance) / sinDistance;
+
+        double x = A * cosLat1 * cos(lon1radians) + B * cosLat2 * cos(lon2radians);
+        double y = A * cosLat1 * sin(lon1radians) + B * cosLat2 * sin(lon2radians);
+        double z = A * sin(lat1radians) + B * sin(lat2radians);
+
+        double midLat = Math.toDegrees(Math.atan2(z, Math.sqrt(x*x + y*y)));
+        double midLon = Math.toDegrees(Math.atan2(y, x));
+
+        return new GHPoint(midLat, midLon);
+    }
+
     @Override
     public boolean isCrossBoundary(double lon1, double lon2) {
         return abs(lon1 - lon2) > 300;
diff --git a/api/src/main/java/com/graphhopper/util/Parameters.java b/api/src/main/java/com/graphhopper/util/Parameters.java
index c927c80ed..5726f8b08 100644
--- a/api/src/main/java/com/graphhopper/util/Parameters.java
+++ b/api/src/main/java/com/graphhopper/util/Parameters.java
@@ -110,6 +110,7 @@ public class Parameters {
          */
         public static final String WAY_POINT_MAX_DISTANCE = "way_point_max_distance";
         public static final String INIT_WAY_POINT_MAX_DISTANCE = ROUTING_INIT_PREFIX + "way_point_max_distance";
+        public static final String ELEVATION_WAY_POINT_MAX_DISTANCE = "elevation_way_point_max_distance";
         /**
          * true or false. If routes at via points should avoid u-turns. (not for CH) See related
          * 'heading' parameter:
diff --git a/config-example.yml b/config-example.yml
index f9447ef83..22663b6eb 100644
--- a/config-example.yml
+++ b/config-example.yml
@@ -92,6 +92,19 @@ graphhopper:
   # graph.elevation.dataaccess: RAM_STORE
 
 
+  # To enable bilinear interpolation when sampling elevation at points (default uses nearest neighbor):
+  # graph.elevation.interpolate: bilinear
+
+
+  # To increase elevation profile resolution, use the following two parameters to tune the extra resolution you need
+  # against the additional storage space used for edge geometries. You should enable bilinear interpolation when using
+  # these features (see #1953 for details).
+  # - first, set the distance (in meters) at which elevation samples should be taken on long edges
+  # graph.elevation.long_edge_sampling_distance: 60
+  # - second, set the elevation tolerance (in meters) to use when simplifying polylines since the default ignores
+  #   elevation and will remove the extra points that long edge sampling added
+  # graph.elevation.way_point_max_distance: 10
+
 
   #### Speed, hybrid and flexible mode ####
 
diff --git a/core/src/main/java/com/graphhopper/GraphHopper.java b/core/src/main/java/com/graphhopper/GraphHopper.java
index b8e7c817a..3f70251e7 100644
--- a/core/src/main/java/com/graphhopper/GraphHopper.java
+++ b/core/src/main/java/com/graphhopper/GraphHopper.java
@@ -86,6 +86,7 @@ public class GraphHopper implements GraphHopperAPI {
     private boolean allowWrites = true;
     private boolean fullyLoaded = false;
     private boolean smoothElevation = false;
+    private double longEdgeSamplingDistance = Double.MAX_VALUE;
     // for routing
     private final RouterConfig routerConfig = new RouterConfig();
     // for index
@@ -328,6 +329,22 @@ public class GraphHopper implements GraphHopperAPI {
         return this;
     }
 
+    /**
+     * Sets the distance distance between elevation samples on long edges
+     */
+    public GraphHopper setLongEdgeSamplingDistance(double longEdgeSamplingDistance) {
+        this.longEdgeSamplingDistance = longEdgeSamplingDistance;
+        return this;
+    }
+
+    /**
+     * Sets the max elevation discrepancy between way points and the simplified polyline in meters
+     */
+    public GraphHopper setElevationWayPointMaxDistance(double elevationWayPointMaxDistance) {
+        this.routerConfig.setElevationWayPointMaxDistance(elevationWayPointMaxDistance);
+        return this;
+    }
+
     public String getGraphHopperLocation() {
         return ghLocation;
     }
@@ -485,9 +502,14 @@ public class GraphHopper implements GraphHopperAPI {
 
         // elevation
         this.smoothElevation = ghConfig.getBool("graph.elevation.smoothing", false);
+        this.longEdgeSamplingDistance = ghConfig.getDouble("graph.elevation.long_edge_sampling_distance", Double.MAX_VALUE);
+        setElevationWayPointMaxDistance(ghConfig.getDouble("graph.elevation.way_point_max_distance", Double.MAX_VALUE));
         ElevationProvider elevationProvider = createElevationProvider(ghConfig);
         setElevationProvider(elevationProvider);
 
+        if (longEdgeSamplingDistance < Double.MAX_VALUE && !elevationProvider.getInterpolate())
+            logger.warn("Long edge sampling enabled, but bilinear interpolation disabled. See #1953");
+
         // optimizable prepare
         minNetworkSize = ghConfig.getInt("prepare.min_network_size", minNetworkSize);
 
@@ -687,7 +709,9 @@ public class GraphHopper implements GraphHopperAPI {
                 setElevationProvider(eleProvider).
                 setWorkerThreads(dataReaderWorkerThreads).
                 setWayPointMaxDistance(dataReaderWayPointMaxDistance).
-                setSmoothElevation(this.smoothElevation);
+                setWayPointElevationMaxDistance(routerConfig.getElevationWayPointMaxDistance()).
+                setSmoothElevation(smoothElevation).
+                setLongEdgeSamplingDistance(longEdgeSamplingDistance);
     }
 
     /**
diff --git a/core/src/main/java/com/graphhopper/reader/DataReader.java b/core/src/main/java/com/graphhopper/reader/DataReader.java
index 70d4282a2..0d7d2f7b2 100644
--- a/core/src/main/java/com/graphhopper/reader/DataReader.java
+++ b/core/src/main/java/com/graphhopper/reader/DataReader.java
@@ -36,8 +36,12 @@ public interface DataReader {
 
     DataReader setWayPointMaxDistance(double wayPointMaxDistance);
 
+    DataReader setWayPointElevationMaxDistance(double elevationWayPointMaxDistance);
+
     DataReader setSmoothElevation(boolean smoothElevation);
 
+    DataReader setLongEdgeSamplingDistance(double longEdgeSamplingDistance);
+
     /**
      * This method triggers reading the underlying data to create a graph
      */
diff --git a/core/src/main/java/com/graphhopper/reader/dem/AbstractElevationProvider.java b/core/src/main/java/com/graphhopper/reader/dem/AbstractElevationProvider.java
index 676eab263..e96a07ef4 100644
--- a/core/src/main/java/com/graphhopper/reader/dem/AbstractElevationProvider.java
+++ b/core/src/main/java/com/graphhopper/reader/dem/AbstractElevationProvider.java
@@ -59,6 +59,11 @@ public abstract class AbstractElevationProvider implements ElevationProvider {
         this.interpolate = interpolate;
     }
 
+    @Override
+    public boolean getInterpolate() {
+        return this.interpolate;
+    }
+
     void setSleep(long sleep) {
         this.sleep = sleep;
     }
diff --git a/core/src/main/java/com/graphhopper/reader/dem/EdgeSampling.java b/core/src/main/java/com/graphhopper/reader/dem/EdgeSampling.java
new file mode 100644
index 000000000..29b206b03
--- /dev/null
+++ b/core/src/main/java/com/graphhopper/reader/dem/EdgeSampling.java
@@ -0,0 +1,77 @@
+/*
+ *  Licensed to GraphHopper GmbH under one or more contributor
+ *  license agreements. See the NOTICE file distributed with this work for
+ *  additional information regarding copyright ownership.
+ *
+ *  GraphHopper GmbH licenses this file to you under the Apache License,
+ *  Version 2.0 (the "License"); you may not use this file except in
+ *  compliance with the License. You may obtain a copy of the License at
+ *
+ *       http://www.apache.org/licenses/LICENSE-2.0
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ */
+package com.graphhopper.reader.dem;
+
+import com.graphhopper.util.DistanceCalc;
+import com.graphhopper.util.DistanceCalcEarth;
+import com.graphhopper.util.PointList;
+import com.graphhopper.util.shapes.GHPoint;
+
+/**
+ * Ensures that elevation is sampled along a point list with no more than maxDistance between samples. Works by adding
+ * points along long edges and fetching elevation at each inserted point.
+ *
+ * For short distances this uses a simple linear approximation to interpolate between points and for longer distances it
+ * uses great circle interpolation.
+ */
+public class EdgeSampling {
+    private static final double GREAT_CIRCLE_SEGMENT_LENGTH = DistanceCalcEarth.METERS_PER_DEGREE / 4;
+
+    private EdgeSampling() {}
+
+    public static PointList sample(long wayOsmId, PointList input, double maxDistance, DistanceCalc distCalc, ElevationProvider elevation) {
+        PointList output = new PointList(input.getSize() * 2, input.is3D());
+        if (input.isEmpty()) return output;
+        int nodes = input.getSize();
+        double lastLat = input.getLat(0), lastLon = input.getLon(0), lastEle = input.getEle(0),
+                thisLat, thisLon, thisEle;
+        for (int i = 0; i < nodes; i++) {
+            thisLat = input.getLat(i);
+            thisLon = input.getLon(i);
+            thisEle = input.getEle(i);
+            if (i > 0) {
+                double segmentLength = distCalc.calcDist3D(lastLat, lastLon, lastEle, thisLat, thisLon, thisEle);
+                int segments = (int) Math.round(segmentLength / maxDistance);
+                // for small distances, we use a simple and fast approximation to interpolate between points
+                // for longer distances (or when crossing international date line) we use great circle interpolation
+                boolean exact = segmentLength > GREAT_CIRCLE_SEGMENT_LENGTH || distCalc.isCrossBoundary(lastLon, thisLon);
+                for (int segment = 1; segment < segments; segment++) {
+                    double ratio = (double) segment / segments;
+                    double lat, lon;
+                    if (exact) {
+                        GHPoint point = distCalc.intermediatePoint(ratio, lastLat, lastLon, thisLat, thisLon);
+                        lat = point.getLat();
+                        lon = point.getLon();
+                    } else {
+                        lat = lastLat + (thisLat - lastLat) * ratio;
+                        lon = lastLon + (thisLon - lastLon) * ratio;
+                    }
+                    double ele = elevation.getEle(lat, lon);
+                    if (!Double.isNaN(ele)) {
+                        output.add(lat, lon, ele);
+                    }
+                }
+            }
+            output.add(thisLat, thisLon, thisEle);
+            lastLat = thisLat;
+            lastLon = thisLon;
+            lastEle = thisEle;
+        }
+        return output;
+    }
+}
diff --git a/core/src/main/java/com/graphhopper/reader/dem/ElevationProvider.java b/core/src/main/java/com/graphhopper/reader/dem/ElevationProvider.java
index 7fc3fc7f0..156eea1f5 100644
--- a/core/src/main/java/com/graphhopper/reader/dem/ElevationProvider.java
+++ b/core/src/main/java/com/graphhopper/reader/dem/ElevationProvider.java
@@ -50,6 +50,11 @@ public interface ElevationProvider {
         @Override
         public void setInterpolate(boolean interpolate) {
         }
+
+        @Override
+        public boolean getInterpolate() {
+            return false;
+        }
     };
 
     /**
@@ -75,6 +80,11 @@ public interface ElevationProvider {
      */
     void setInterpolate(boolean interpolate);
 
+    /**
+     * Returns true if bilinear interpolation is enabled.
+     */
+    boolean getInterpolate();
+
     /**
      * Release resources.
      */
diff --git a/core/src/main/java/com/graphhopper/reader/dem/MultiSourceElevationProvider.java b/core/src/main/java/com/graphhopper/reader/dem/MultiSourceElevationProvider.java
index 451a7d4ff..fc6937dd3 100644
--- a/core/src/main/java/com/graphhopper/reader/dem/MultiSourceElevationProvider.java
+++ b/core/src/main/java/com/graphhopper/reader/dem/MultiSourceElevationProvider.java
@@ -82,6 +82,11 @@ public class MultiSourceElevationProvider implements ElevationProvider {
         globalProvider.setInterpolate(interpolate);
     }
 
+    @Override
+    public boolean getInterpolate() {
+        return srtmProvider.getInterpolate() && globalProvider.getInterpolate();
+    }
+
     @Override
     public void release() {
         srtmProvider.release();
diff --git a/core/src/main/java/com/graphhopper/routing/Router.java b/core/src/main/java/com/graphhopper/routing/Router.java
index 007344585..01a6c8297 100644
--- a/core/src/main/java/com/graphhopper/routing/Router.java
+++ b/core/src/main/java/com/graphhopper/routing/Router.java
@@ -279,8 +279,11 @@ public class Router {
         boolean enableInstructions = request.getHints().getBool(Parameters.Routing.INSTRUCTIONS, encodingManager.isEnableInstructions());
         boolean calcPoints = request.getHints().getBool(Parameters.Routing.CALC_POINTS, routerConfig.isCalcPoints());
         double wayPointMaxDistance = request.getHints().getDouble(Parameters.Routing.WAY_POINT_MAX_DISTANCE, 1d);
+        double elevationWayPointMaxDistance = request.getHints().getDouble(ELEVATION_WAY_POINT_MAX_DISTANCE, routerConfig.getElevationWayPointMaxDistance());
 
-        DouglasPeucker peucker = new DouglasPeucker().setMaxDistance(wayPointMaxDistance);
+        DouglasPeucker peucker = new DouglasPeucker().
+                setMaxDistance(wayPointMaxDistance).
+                setElevationMaxDistance(elevationWayPointMaxDistance);
         PathMerger pathMerger = new PathMerger(graph, weighting).
                 setCalcPoints(calcPoints).
                 setDouglasPeucker(peucker).
diff --git a/core/src/main/java/com/graphhopper/routing/RouterConfig.java b/core/src/main/java/com/graphhopper/routing/RouterConfig.java
index aeb73d933..de7e4ca7c 100644
--- a/core/src/main/java/com/graphhopper/routing/RouterConfig.java
+++ b/core/src/main/java/com/graphhopper/routing/RouterConfig.java
@@ -29,6 +29,7 @@ public class RouterConfig {
     private boolean lmDisablingAllowed = true;
     private boolean calcPoints = true;
     private boolean simplifyResponse = true;
+    private double elevationWayPointMaxDistance = Double.MAX_VALUE;
     private int activeLandmarkCount = 8;
 
     public int getMaxVisitedNodes() {
@@ -105,4 +106,12 @@ public class RouterConfig {
     public void setActiveLandmarkCount(int activeLandmarkCount) {
         this.activeLandmarkCount = activeLandmarkCount;
     }
+
+    public double getElevationWayPointMaxDistance() {
+        return elevationWayPointMaxDistance;
+    }
+
+    public void setElevationWayPointMaxDistance(double elevationWayPointMaxDistance) {
+        this.elevationWayPointMaxDistance = elevationWayPointMaxDistance;
+    }
 }
diff --git a/core/src/main/java/com/graphhopper/util/DistanceCalcEuclidean.java b/core/src/main/java/com/graphhopper/util/DistanceCalcEuclidean.java
index 276ab1f02..5e6acb2a7 100644
--- a/core/src/main/java/com/graphhopper/util/DistanceCalcEuclidean.java
+++ b/core/src/main/java/com/graphhopper/util/DistanceCalcEuclidean.java
@@ -92,6 +92,15 @@ public class DistanceCalcEuclidean extends DistanceCalcEarth {
         throw new UnsupportedOperationException("Not supported for the 2D Euclidean space");
     }
 
+    @Override
+    public GHPoint intermediatePoint(double f, double lat1, double lon1, double lat2, double lon2) {
+        double delatLat = lat2 - lat1;
+        double deltaLon = lon2 - lon1;
+        double midLat = lat1 + delatLat * f;
+        double midLon = lon1 + deltaLon * f;
+        return new GHPoint(midLat, midLon);
+    }
+
     @Override
     public boolean isCrossBoundary(double lon1, double lon2) {
         throw new UnsupportedOperationException("Not supported for the 2D Euclidean space");
@@ -118,15 +127,16 @@ public class DistanceCalcEuclidean extends DistanceCalcEarth {
 
         double norm = dx * dx + dy * dy + dz * dz;
         double factor = ((rx - ax) * dx + (ry - ay) * dy + (rz - az) * dz) / norm;
+        if (Double.isNaN(factor)) factor = 0;
 
         // x,y,z is projection of r onto segment a-b
         double cx = ax + factor * dx;
         double cy = ay + factor * dy;
         double cz = az + factor * dz;
 
-        double rdx = cx - ax;
-        double rdy = cy - ay;
-        double rdz = cz - az;
+        double rdx = cx - rx;
+        double rdy = cy - ry;
+        double rdz = cz - rz;
 
         return rdx * rdx + rdy * rdy + rdz * rdz;
     }
diff --git a/core/src/main/java/com/graphhopper/util/DouglasPeucker.java b/core/src/main/java/com/graphhopper/util/DouglasPeucker.java
index c7acdfdc3..90e889676 100644
--- a/core/src/main/java/com/graphhopper/util/DouglasPeucker.java
+++ b/core/src/main/java/com/graphhopper/util/DouglasPeucker.java
@@ -28,6 +28,8 @@ package com.graphhopper.util;
  */
 public class DouglasPeucker {
     private double normedMaxDist;
+    private double elevationMaxDistance;
+    private double maxDistance;
     private DistanceCalc calc;
     private boolean approx;
 
@@ -35,6 +37,8 @@ public class DouglasPeucker {
         setApproximation(true);
         // 1m
         setMaxDistance(1);
+        // elevation ignored by default
+        setElevationMaxDistance(Double.MAX_VALUE);
     }
 
     public void setApproximation(boolean a) {
@@ -50,6 +54,15 @@ public class DouglasPeucker {
      */
     public DouglasPeucker setMaxDistance(double dist) {
         this.normedMaxDist = calc.calcNormalizedDist(dist);
+        this.maxDistance = dist;
+        return this;
+    }
+
+    /**
+     * maximum elevation distance of discrepancy (from the normal way) in meters
+     */
+    public DouglasPeucker setElevationMaxDistance(double dist) {
+        this.elevationMaxDistance = dist;
         return this;
     }
 
@@ -108,17 +121,26 @@ public class DouglasPeucker {
         }
         int indexWithMaxDist = -1;
         double maxDist = -1;
+        double elevationFactor = maxDistance / elevationMaxDistance;
         double firstLat = points.getLatitude(fromIndex);
         double firstLon = points.getLongitude(fromIndex);
+        double firstEle = points.getElevation(fromIndex);
         double lastLat = points.getLatitude(lastIndex);
         double lastLon = points.getLongitude(lastIndex);
+        double lastEle = points.getElevation(lastIndex);
         for (int i = fromIndex + 1; i < lastIndex; i++) {
             double lat = points.getLatitude(i);
             if (Double.isNaN(lat)) {
                 continue;
             }
             double lon = points.getLongitude(i);
-            double dist = calc.calcNormalizedEdgeDistance(lat, lon, firstLat, firstLon, lastLat, lastLon);
+            double ele = points.getElevation(i);
+            double dist = (points.is3D() && elevationMaxDistance < Double.MAX_VALUE && !Double.isNaN(firstEle) && !Double.isNaN(lastEle) && !Double.isNaN(ele))
+                    ? calc.calcNormalizedEdgeDistance3D(
+                            lat, lon, ele * elevationFactor,
+                            firstLat, firstLon, firstEle * elevationFactor,
+                            lastLat, lastLon, lastEle * elevationFactor)
+                    : calc.calcNormalizedEdgeDistance(lat, lon, firstLat, firstLon, lastLat, lastLon);
             if (maxDist < dist) {
                 indexWithMaxDist = i;
                 maxDist = dist;
diff --git a/reader-gtfs/src/main/java/com/graphhopper/reader/gtfs/GraphHopperGtfs.java b/reader-gtfs/src/main/java/com/graphhopper/reader/gtfs/GraphHopperGtfs.java
index 5d5abc1f3..a97fc89b3 100644
--- a/reader-gtfs/src/main/java/com/graphhopper/reader/gtfs/GraphHopperGtfs.java
+++ b/reader-gtfs/src/main/java/com/graphhopper/reader/gtfs/GraphHopperGtfs.java
@@ -92,11 +92,21 @@ public class GraphHopperGtfs extends GraphHopperOSM {
                     return this;
                 }
 
+                @Override
+                public DataReader setWayPointElevationMaxDistance(double elevationWayPointMaxDistance) {
+                    return this;
+                }
+
                 @Override
                 public DataReader setSmoothElevation(boolean smoothElevation) {
                     return this;
                 }
 
+                @Override
+                public DataReader setLongEdgeSamplingDistance(double longEdgeSamplingDistance) {
+                    return this;
+                }
+
                 @Override
                 public void readGraph() {
 
diff --git a/reader-osm/src/main/java/com/graphhopper/reader/osm/OSMReader.java b/reader-osm/src/main/java/com/graphhopper/reader/osm/OSMReader.java
index 932a0701b..3859865ff 100644
--- a/reader-osm/src/main/java/com/graphhopper/reader/osm/OSMReader.java
+++ b/reader-osm/src/main/java/com/graphhopper/reader/osm/OSMReader.java
@@ -21,6 +21,7 @@ import com.carrotsearch.hppc.*;
 import com.graphhopper.coll.LongIntMap;
 import com.graphhopper.coll.*;
 import com.graphhopper.reader.*;
+import com.graphhopper.reader.dem.EdgeSampling;
 import com.graphhopper.reader.dem.ElevationProvider;
 import com.graphhopper.reader.dem.GraphElevationSmoothing;
 import com.graphhopper.routing.ev.BooleanEncodedValue;
@@ -77,6 +78,7 @@ public class OSMReader implements DataReader, TurnCostParser.ExternalInternalMap
     private final DistanceCalc distCalc = Helper.DIST_EARTH;
     private final DouglasPeucker simplifyAlgo = new DouglasPeucker();
     private boolean smoothElevation = false;
+    private double longEdgeSamplingDistance = 0;
     private final boolean exitOnlyPillarNodeException = true;
     private final Map<String, EdgeExplorer> outExplorerMap = new HashMap<>();
     private final Map<String, EdgeExplorer> inExplorerMap = new HashMap<>();
@@ -670,6 +672,10 @@ public class OSMReader implements DataReader, TurnCostParser.ExternalInternalMap
         if (this.smoothElevation)
             pointList = GraphElevationSmoothing.smoothElevation(pointList);
 
+        // sample points along long edges
+        if (this.longEdgeSamplingDistance < Double.MAX_VALUE && pointList.is3D())
+            pointList = EdgeSampling.sample(wayOsmId, pointList, longEdgeSamplingDistance, distCalc, eleProvider);
+
         double towerNodeDistance = pointList.calcDistance(distCalc);
 
         if (towerNodeDistance < 0.001) {
@@ -896,12 +902,24 @@ public class OSMReader implements DataReader, TurnCostParser.ExternalInternalMap
         return this;
     }
 
+    @Override
+    public DataReader setWayPointElevationMaxDistance(double elevationWayPointMaxDistance) {
+        simplifyAlgo.setElevationMaxDistance(elevationWayPointMaxDistance);
+        return this;
+    }
+
     @Override
     public DataReader setSmoothElevation(boolean smoothElevation) {
         this.smoothElevation = smoothElevation;
         return this;
     }
 
+    @Override
+    public DataReader setLongEdgeSamplingDistance(double longEdgeSamplingDistance) {
+        this.longEdgeSamplingDistance = longEdgeSamplingDistance;
+        return this;
+    }
+
     @Override
     public OSMReader setWorkerThreads(int numOfWorkers) {
         this.workerThreads = numOfWorkers;
diff --git a/web-bundle/src/main/java/com/graphhopper/resources/RouteResource.java b/web-bundle/src/main/java/com/graphhopper/resources/RouteResource.java
index f46327c7f..a100459f9 100644
--- a/web-bundle/src/main/java/com/graphhopper/resources/RouteResource.java
+++ b/web-bundle/src/main/java/com/graphhopper/resources/RouteResource.java
@@ -73,6 +73,7 @@ public class RouteResource {
             @Context HttpServletRequest httpReq,
             @Context UriInfo uriInfo,
             @QueryParam(WAY_POINT_MAX_DISTANCE) @DefaultValue("1") double minPathPrecision,
+            @QueryParam(ELEVATION_WAY_POINT_MAX_DISTANCE) Double minPathElevationPrecision,
             @QueryParam("point") @NotNull List<GHPointParam> pointParams,
             @QueryParam("type") @DefaultValue("json") String type,
             @QueryParam(INSTRUCTIONS) @DefaultValue("true") boolean instructions,
@@ -122,6 +123,10 @@ public class RouteResource {
                 putObject(INSTRUCTIONS, instructions).
                 putObject(WAY_POINT_MAX_DISTANCE, minPathPrecision);
 
+        if (minPathElevationPrecision != null) {
+            request.getHints().putObject(ELEVATION_WAY_POINT_MAX_DISTANCE, minPathElevationPrecision);
+        }
+
         GHResponse ghResponse = graphHopper.route(request);
 
         long took = sw.stop().getNanos() / 1_000_000;


Test patch:
diff --git a/api/src/test/java/com/graphhopper/util/DistanceCalcEarthTest.java b/api/src/test/java/com/graphhopper/util/DistanceCalcEarthTest.java
index 921ad97d1..7acf1596a 100644
--- a/api/src/test/java/com/graphhopper/util/DistanceCalcEarthTest.java
+++ b/api/src/test/java/com/graphhopper/util/DistanceCalcEarthTest.java
@@ -142,6 +142,46 @@ public class DistanceCalcEarthTest {
         assertEquals(10, calc.calcDenormalizedDist(dist), 1e-4);
     }
 
+    @Test
+    public void testEdgeDistanceStartEndSame() {
+        DistanceCalc calc = new DistancePlaneProjection();
+        // just change elevation
+        double dist = calc.calcNormalizedEdgeDistance3D(0, 0, 10,
+                0, 0, 0,
+                0, 0, 0);
+        assertEquals(10, calc.calcDenormalizedDist(dist), 1e-4);
+        // just change lat
+        dist = calc.calcNormalizedEdgeDistance3D(1, 0, 0,
+                0, 0, 0,
+                0, 0, 0);
+        assertEquals(DistanceCalcEarth.METERS_PER_DEGREE, calc.calcDenormalizedDist(dist), 1e-4);
+        // just change lon
+        dist = calc.calcNormalizedEdgeDistance3D(0, 1, 0,
+                0, 0, 0,
+                0, 0, 0);
+        assertEquals(DistanceCalcEarth.METERS_PER_DEGREE, calc.calcDenormalizedDist(dist), 1e-4);
+    }
+
+    @Test
+    public void testEdgeDistanceStartEndDifferentElevation() {
+        DistanceCalc calc = new DistancePlaneProjection();
+        // just change elevation
+        double dist = calc.calcNormalizedEdgeDistance3D(0, 0, 10,
+                0, 0, 0,
+                0, 0, 1);
+        assertEquals(0, calc.calcDenormalizedDist(dist), 1e-4);
+        // just change lat
+        dist = calc.calcNormalizedEdgeDistance3D(1, 0, 0,
+                0, 0, 0,
+                0, 0, 1);
+        assertEquals(DistanceCalcEarth.METERS_PER_DEGREE, calc.calcDenormalizedDist(dist), 1e-4);
+        // just change lon
+        dist = calc.calcNormalizedEdgeDistance3D(0, 1, 0,
+                0, 0, 0,
+                0, 0, 1);
+        assertEquals(DistanceCalcEarth.METERS_PER_DEGREE, calc.calcDenormalizedDist(dist), 1e-4);
+    }
+
     @Test
     public void testValidEdgeDistance() {
         assertTrue(dc.validEdgeDistance(49.94241, 11.544356, 49.937964, 11.541824, 49.942272, 11.555643));
@@ -257,4 +297,38 @@ public class DistanceCalcEarthTest {
                 0, 0, Double.NaN
         ), 1e-6);
     }
+
+    @Test
+    public void testIntermediatePoint() {
+        DistanceCalc distCalc = new DistanceCalcEarth();
+        GHPoint point = distCalc.intermediatePoint(0, 0, 0, 0, 0);
+        assertEquals(0, point.getLat(), 1e-5);
+        assertEquals(0, point.getLon(), 1e-5);
+
+        point = distCalc.intermediatePoint(0.5, 0, 0, 10, 0);
+        assertEquals(5, point.getLat(), 1e-5);
+        assertEquals(0, point.getLon(), 1e-5);
+
+        point = distCalc.intermediatePoint(0.5, 0, 0, 0, 10);
+        assertEquals(0, point.getLat(), 1e-5);
+        assertEquals(5, point.getLon(), 1e-5);
+
+        // cross international date line going west
+        point = distCalc.intermediatePoint(0.5, 45, -179, 45, 177);
+        assertEquals(45, point.getLat(), 1);
+        assertEquals(179, point.getLon(), 1e-5);
+
+        // cross international date line going east
+        point = distCalc.intermediatePoint(0.5, 45, 179, 45, -177);
+        assertEquals(45, point.getLat(), 1);
+        assertEquals(-179, point.getLon(), 1e-5);
+
+        // cross north pole
+        point = distCalc.intermediatePoint(0.25, 45, -90, 45, 90);
+        assertEquals(67.5, point.getLat(), 1e-1);
+        assertEquals(-90, point.getLon(), 1e-5);
+        point = distCalc.intermediatePoint(0.75, 45, -90, 45, 90);
+        assertEquals(67.5, point.getLat(), 1e-1);
+        assertEquals(90, point.getLon(), 1e-5);
+    }
 }
diff --git a/core/src/test/java/com/graphhopper/reader/dem/EdgeSamplingTest.java b/core/src/test/java/com/graphhopper/reader/dem/EdgeSamplingTest.java
new file mode 100644
index 000000000..75c8131ac
--- /dev/null
+++ b/core/src/test/java/com/graphhopper/reader/dem/EdgeSamplingTest.java
@@ -0,0 +1,137 @@
+package com.graphhopper.reader.dem;
+
+import com.graphhopper.util.DistanceCalcEarth;
+import com.graphhopper.util.PointList;
+import org.junit.Test;
+
+import static org.junit.Assert.assertEquals;
+
+public class EdgeSamplingTest {
+    private final ElevationProvider elevation = new AbstractElevationProvider("") {
+        @Override
+        public double getEle(double lat, double lon) {
+            return 10;
+        }
+
+        @Override
+        public void release() {}
+        @Override
+        String getFileName(double lat, double lon) { return ""; }
+        @Override
+        String getDownloadURL(double lat, double lon) { return ""; }
+    };
+
+    private double round(double d) {
+        return Math.round(d * 1000) / 1000.0;
+    }
+
+    private PointList round(PointList list) {
+        for (int i = 0; i < list.size(); i++) {
+            list.set(i, round(list.getLat(i)), round(list.getLon(i)), list.getEle(i));
+        }
+        return list;
+    }
+
+    @Test
+    public void doesNotAddExtraPointBelowThreshold() {
+        PointList in = new PointList(2, true);
+        in.add(0, 0, 0);
+        in.add(1.4, 0, 0);
+
+        PointList out = EdgeSampling.sample(
+                1,
+                in,
+                DistanceCalcEarth.METERS_PER_DEGREE,
+                new DistanceCalcEarth(),
+                elevation
+        );
+
+        assertEquals("(0.0,0.0,0.0), (1.4,0.0,0.0)", round(out).toString());
+    }
+
+    @Test
+    public void addsExtraPointAboveThreshold() {
+        PointList in = new PointList(2, true);
+        in.add(0, 0, 0);
+        in.add(0.8, 0, 0);
+
+        PointList out = EdgeSampling.sample(
+                1,
+                in,
+                DistanceCalcEarth.METERS_PER_DEGREE / 2,
+                new DistanceCalcEarth(),
+                elevation
+        );
+
+        assertEquals("(0.0,0.0,0.0), (0.4,0.0,10.0), (0.8,0.0,0.0)", round(out).toString());
+    }
+
+    @Test
+    public void addsExtraPointBelowSecondThreshold() {
+        PointList in = new PointList(2, true);
+        in.add(0, 0, 0);
+        in.add(0.8, 0, 0);
+
+        PointList out = EdgeSampling.sample(
+                1,
+                in,
+                DistanceCalcEarth.METERS_PER_DEGREE / 3,
+                new DistanceCalcEarth(),
+                elevation
+        );
+
+        assertEquals("(0.0,0.0,0.0), (0.4,0.0,10.0), (0.8,0.0,0.0)", round(out).toString());
+    }
+
+    @Test
+    public void addsTwoPointsAboveThreshold() {
+        PointList in = new PointList(2, true);
+        in.add(0, 0, 0);
+        in.add(0.75, 0, 0);
+
+        PointList out = EdgeSampling.sample(
+                1,
+                in,
+                DistanceCalcEarth.METERS_PER_DEGREE / 4,
+                new DistanceCalcEarth(),
+                elevation
+        );
+
+        assertEquals("(0.0,0.0,0.0), (0.25,0.0,10.0), (0.5,0.0,10.0), (0.75,0.0,0.0)", round(out).toString());
+    }
+
+    @Test
+    public void doesntAddPointsCrossingInternationalDateLine() {
+        PointList in = new PointList(2, true);
+        in.add(0, -178.5, 0);
+        in.add(0.0, 178.5, 0);
+
+        PointList out = EdgeSampling.sample(
+                1,
+                in,
+                DistanceCalcEarth.METERS_PER_DEGREE,
+                new DistanceCalcEarth(),
+                elevation
+        );
+
+        assertEquals("(0.0,-178.5,0.0), (0.0,-179.5,10.0), (0.0,179.5,10.0), (0.0,178.5,0.0)", round(out).toString());
+    }
+
+    @Test
+    public void usesGreatCircleInterpolationOnLongPaths() {
+        PointList in = new PointList(2, true);
+        in.add(88.5, -90, 0);
+        in.add(88.5, 90, 0);
+
+        PointList out = EdgeSampling.sample(
+                1,
+                in,
+                DistanceCalcEarth.METERS_PER_DEGREE,
+                new DistanceCalcEarth(),
+                elevation
+        );
+
+        assertEquals("(88.5,-90.0,0.0), (89.5,-90.0,10.0), (89.5,90.0,10.0), (88.5,90.0,0.0)", round(out).toString());
+    }
+
+}
\ No newline at end of file
diff --git a/core/src/test/java/com/graphhopper/util/DistanceCalcEuclideanTest.java b/core/src/test/java/com/graphhopper/util/DistanceCalcEuclideanTest.java
index fca9d7cbe..f6502528b 100644
--- a/core/src/test/java/com/graphhopper/util/DistanceCalcEuclideanTest.java
+++ b/core/src/test/java/com/graphhopper/util/DistanceCalcEuclideanTest.java
@@ -40,6 +40,13 @@ public class DistanceCalcEuclideanTest {
         Assert.assertEquals(50, distance, 0);
     }
 
+    @Test
+    public void testCalcNormalizedEdgeDistance3dStartEndSame() {
+        DistanceCalcEuclidean distanceCalc = new DistanceCalcEuclidean();
+        double distance = distanceCalc.calcNormalizedEdgeDistance3D(0, 3, 4, 0, 0, 0, 0, 0, 0);
+        Assert.assertEquals(25, distance, 0);
+    }
+
     @Test
     public void testValidEdgeDistance() {
         DistanceCalcEuclidean distanceCalc = new DistanceCalcEuclidean();
diff --git a/core/src/test/java/com/graphhopper/util/DouglasPeuckerTest.java b/core/src/test/java/com/graphhopper/util/DouglasPeuckerTest.java
index cae31b382..9f0feabcc 100644
--- a/core/src/test/java/com/graphhopper/util/DouglasPeuckerTest.java
+++ b/core/src/test/java/com/graphhopper/util/DouglasPeuckerTest.java
@@ -137,4 +137,69 @@ public class DouglasPeuckerTest {
         }
         assertEquals(expected, given);
     }
+
+    @Test
+    public void test3dPathSimplify() {
+        PointList pointList = new PointList(5, true);
+        pointList.add(0, 0, 0);
+        pointList.add(0.01, 0, 10); // can be removed
+        pointList.add(0.02, 0, 20); // can be removed
+        pointList.add(0.03, 0, 30); // can't be removed
+        pointList.add(0.04, 0, 50);
+        new DouglasPeucker().setMaxDistance(1).setElevationMaxDistance(1).simplify(pointList);
+        assertEquals("(0.0,0.0,0.0), (0.03,0.0,30.0), (0.04,0.0,50.0)", pointList.toString());
+    }
+
+    @Test
+    public void test3dPathSimplifyElevationDisabled() {
+        PointList pointList = new PointList(5, true);
+        pointList.add(0, 0, 0);
+        pointList.add(0.03, 0, 30); // would be kept, if we cared about elevation
+        pointList.add(0.04, 0, 50);
+        new DouglasPeucker().setMaxDistance(1).setElevationMaxDistance(Double.MAX_VALUE).simplify(pointList);
+        assertEquals("(0.0,0.0,0.0), (0.04,0.0,50.0)", pointList.toString());
+    }
+
+    @Test
+    public void test3dPathSimplifyElevationMaxDistFive() {
+        PointList pointList = new PointList(5, true);
+        pointList.add(0, 0, 0);
+        pointList.add(0.01, 0, 14); // <5m from straight line (10), remove
+        pointList.add(0.02, 0, 20); // on straight line, remove
+        pointList.add(0.03, 0, 30); // >5m from straight line, keep
+        pointList.add(0.04, 0, 50);
+        new DouglasPeucker().setMaxDistance(1).setElevationMaxDistance(5).simplify(pointList);
+        assertEquals("(0.0,0.0,0.0), (0.03,0.0,30.0), (0.04,0.0,50.0)", pointList.toString());
+    }
+
+    @Test
+    public void test3dPathSimplifyWithMissingElevation() {
+        PointList pointList = new PointList(5, true);
+        pointList.add(0, 0, 0);
+        pointList.add(0, 0.5, Double.NaN); // on straight line in 2d space, ignore elevation
+        pointList.add(0, 1, 14); // <5m from straight line (10), remove
+        pointList.add(1, 1, 20);
+        new DouglasPeucker().setMaxDistance(1).setElevationMaxDistance(1).simplify(pointList);
+        assertEquals("(0.0,0.0,0.0), (0.0,1.0,14.0), (1.0,1.0,20.0)", pointList.toString());
+    }
+
+    @Test
+    public void test3dSimplifyStartEndSame() {
+        PointList pointList = new PointList(3, true);
+        pointList.add(0, 0, 0);
+        pointList.add(0.03, 0, 30);
+        pointList.add(0, 0, 0);
+        new DouglasPeucker().setMaxDistance(1).setElevationMaxDistance(1).simplify(pointList);
+        assertEquals("(0.0,0.0,0.0), (0.03,0.0,30.0), (0.0,0.0,0.0)", pointList.toString());
+    }
+
+    @Test
+    public void test2dSimplifyStartEndSame() {
+        PointList pointList = new PointList(3, false);
+        pointList.add(0, 0);
+        pointList.add(0.03, 0);
+        pointList.add(0, 0);
+        new DouglasPeucker().setMaxDistance(1).setElevationMaxDistance(1).simplify(pointList);
+        assertEquals("(0.0,0.0), (0.03,0.0), (0.0,0.0)", pointList.toString());
+    }
 }
diff --git a/reader-osm/src/test/java/com/graphhopper/GraphHopperTest.java b/reader-osm/src/test/java/com/graphhopper/GraphHopperTest.java
index 94850535d..aff925ae7 100644
--- a/reader-osm/src/test/java/com/graphhopper/GraphHopperTest.java
+++ b/reader-osm/src/test/java/com/graphhopper/GraphHopperTest.java
@@ -21,6 +21,7 @@ import com.graphhopper.config.CHProfile;
 import com.graphhopper.config.LMProfile;
 import com.graphhopper.config.Profile;
 import com.graphhopper.reader.ReaderWay;
+import com.graphhopper.reader.dem.ElevationProvider;
 import com.graphhopper.reader.dem.SRTMProvider;
 import com.graphhopper.reader.dem.SkadiProvider;
 import com.graphhopper.reader.osm.GraphHopperOSM;
@@ -991,6 +992,50 @@ public class GraphHopperTest {
         assertEquals(31.32, pointList.getEle(5), .1);
     }
 
+    @Test
+    public void testSRTMWithLongEdgeSampling() {
+        final String profile = "profile";
+        final String vehicle = "foot";
+        final String weighting = "shortest";
+
+        GraphHopper hopper = createGraphHopper(vehicle).
+                setOSMFile(MONACO).
+                setStoreOnFlush(true).
+                setElevationWayPointMaxDistance(1).
+                setProfiles(new Profile("profile").setVehicle(vehicle).setWeighting(weighting)).
+                setLongEdgeSamplingDistance(30);
+
+        ElevationProvider elevationProvider = new SRTMProvider(DIR);
+        elevationProvider.setInterpolate(true);
+        hopper.setElevationProvider(elevationProvider);
+        hopper.importOrLoad();
+
+        GHResponse rsp = hopper.route(new GHRequest(43.730729, 7.421288, 43.727697, 7.419199).
+                setAlgorithm(ASTAR).setProfile(profile));
+
+        ResponsePath arsp = rsp.getBest();
+        assertEquals(1570.4, arsp.getDistance(), .1);
+        assertEquals(60, arsp.getPoints().getSize());
+        assertTrue(arsp.getPoints().is3D());
+
+        InstructionList il = arsp.getInstructions();
+        assertEquals(12, il.size());
+        assertTrue(il.get(0).getPoints().is3D());
+
+        String str = arsp.getPoints().toString();
+
+        assertEquals(23.8, arsp.getAscend(), 1e-1);
+        assertEquals(67.4, arsp.getDescend(), 1e-1);
+
+        assertEquals(60, arsp.getPoints().size());
+        assertEquals(new GHPoint3D(43.73068455771767, 7.421283689825812, 55.82900047302246), arsp.getPoints().get(0));
+        assertEquals(new GHPoint3D(43.727679637988224, 7.419198521975086, 12.274499893188477), arsp.getPoints().get(arsp.getPoints().size() - 1));
+
+        assertEquals(55.83, arsp.getPoints().get(0).getElevation(), 1e-2);
+        assertEquals(57.78, arsp.getPoints().get(1).getElevation(), 1e-2);
+        assertEquals(52.43, arsp.getPoints().get(10).getElevation(), 1e-2);
+    }
+
     @Ignore
     public void testSkadiElevationProvider() {
         final String profile = "profile";