Trips incorrectly show elevation as below sea level

[grant-humphries]

I've come across a few trips that return legs or steps as having an elevation below sea level where it seems unlikely that is true. For the two trips below I've checked those stretches against other sources and they show them as being ~40m higher. These we're initially found on TriMet's production instance (v0.10), but I set up a local instance of v0.20 on my PC and the issue persisted.

See examples in the last leg of this trip, and the latter half of this one.

[abyrd]

We looked into this with @demory at our last routing meeting. It appears that this is due to the vertical datum transformation that we apply to the elevation data. It's not clear why this is happening now but didn't happen in the past. The main thing that changed is that the NED data download service has changed so we're using bulk NED data that we copied to S3 once upon a time. Perhaps the vertical datum of that data is different.

[abyrd]

I'm moving some comments over here from Basecamp for the wider OTP community. In #826 (comment) @grant-humphries is converting some NED-compatible LiDAR data of bridge surfaces to WGS84 vertical coordinates and getting negative numbers, about -4 for the Hawthorne Bridge deck and -12 for the lower deck of the Steel Bridge which is right over the water surface, in the river channel.

The WGS84 system has both an ellipsoid (smooth mathematical surface) and a geoid (irregular estimate of mean sea level). WGS84 height coordinates are relative to the ellipsoid, and the geoid can also be expressed in WGS84 height coordinates relative to the ellipsoid. The geoid varies above and below the ellipsoid by design.

Maps of geoid height show that it's 20-30 meters below the ellipsoid in most of the US:

Which is to say that height above sea level is height relative to the WS84 ellipsoid plus 20 something meters, and negative WGS84 heights are completely normal at low elevations above sea level like Portland (15m).

You can also verify this with the geoid calculator at http://www.unavco.org/software/geodetic-utilities/geoid-height-calculator/geoid-height-calculator.html which shows that at the center of Portland, the geoid is about 22 meters below the ellipsoid, and Portland is about 15m above sea level, so you'd be at -7 meters WGS84 altitude standing on the ground. This places the Hawthorne Bridge deck 3 meters above the typical ground elevation of Portland, and the Steel Bridge lower deck 5 meters below, which makes sense to me.

[abyrd]

Long story short, if we report WGS84 elevations they're going to be negative in coastal areas. We should calculate the height of the ellipsoid above the geoid at the center of the location where OTP is deployed and add that to all reported elevations, and/or only ever use elevations in a relative way.

[abyrd]

I have found an implementation of the ellipsoid -> geoid transform in Geotools. It's in the gt-referencing3D module. http://docs.geotools.org/latest/javadocs/org/geotools/referencing/operation/transform/EarthGravitationalModel.html

We should be able to easily transform coordinates all over the world using this.

I am pushing this ticket back to OTP 1.1 since it doesn't appear to be a bug, just an inconvenient discrepancy between GPS height and height above sea level. @demory @fpurcell let's track this issue here from now on because I believe it's of wider interest to the OTP community.

[barbeau]

I have found an implementation of the ellipsoid -> geoid transform in Geotools

@abyrd Thanks for sharing this! It's something I've been looking for.

[abyrd]

After a discussion in the weekly Conveyal OTP check-in we have decided to transform the elevation values at graph build time. We'll use Geotools the EarthGravitationalModel class on every single point, measure speed, if that's too slow we'll just get one value for the center of the region at the beginning of the elevation graph builder module run. This means that all output will be in meters above sea level, not relative to the ellipsoid. Could this cause problems?

[abyrd]

Well rats. It appears that OSM data is supposed to be in meters above the geoid (mean sea level approximation), not relative to the ellipsoid. We set all the bridge heights in Portland relative to the ellipsoid (i.e. some are negative).

From the OSM wiki (http://wiki.openstreetmap.org/wiki/Key:ele):
Elevation (height above sea level) of a point in metres. This is mainly intended for mountain peaks but could also be used for elevation of airport runways and many other objects. For OpenStreetMap, please use the elevation above mean sea level defined by the World Geodetic System, revision WGS 84 (height above WGS84/EGM96 geoid, also known as "corrected" WGS84 height). National "above sea level" systems usually do not differ much from this value (< 1m). This is not the height above the WGS84 ellipsoid (see Geoid) which is shown by some GPS devices.

Considering the minor differences between different local "above sea level" systems, I am wondering why we're converting anything at all. Aren't NED, OSM, LIDAR data etc. all in meters above sea level?

[barbeau]

This means that all output will be in meters above sea level, not relative to the ellipsoid. Could this cause problems?

I think the main downside to this is more confusion by client code, etc., especially when comparing real-time GPS output to the trip plan. How big of a penalty would it be to store both WGS84 ellipsoid and geoid heights?

This would help eliminate confusion and also prevent clients from needing to calculate both themselves.

[abyrd]

@barbeau for all practical purposes, the difference between the two is a constant in a single city. We could store that conversion constant and output it in the response.

[barbeau]

We could store that conversion constant and output it in the response

That works, thanks!

[andreyz]

for all practical purposes, the difference between the two is a constant in a single city

What about instances that are used on a national level? How large should the service area can be for the error to acceptable?

[demory]

@abyrd @barbeau I'm finally getting around to implementing this now, at least the city-specific fix discussed above. Andrew if I understand your proposal correctly, there would be some mechanism for specifying the difference as a constant for a given deployment, and that constant would be included with the response along with the (unadjusted) WGS84 elevations. This means work would be required on the client side to compute and display the "corrected" elevations.

What if the offset constant, when supplied (via router-config.json, presumably), were automatically applied to the elevations returned in the API response? This would be a new configuration parameter so existing deployments would not be affected (they would continue to return the unadjusted heights). I think TriMet is looking for something that could be applied on the server side without a need to update their client(s).

[barbeau]

As long as we're able to derive both altitudes client side, and it's clearly documented what the fields contain, that works for me.

[abyrd]

@demory I think the best approach is in a comment above: "We have decided to transform the elevation values at graph build time. We'll use Geotools the EarthGravitationalModel class on every single point, measure speed, if that's too slow we'll just get one value for the center of the region at the beginning of the elevation graph builder module run."

If you don't want to mess with applying the EarthGravitationalModel class to every point in the elevation builder, just grab a single adjustment constant by applying EarthGravitationalModel to a single point somewhere in the graph (center of the region). That can be added to every number when setting the elevations, and then saved in the graph.

To meet @barbeau's need to have both altitudes client side, we'd need to include the adjustment constant in the responses.

[abyrd]

I think the best way to avoid disrupting existing deployments is to make application of the adjustment constant to street elevations optional, a boolean configuration parameter in the graph build config. Or actually it might be better to do it in the router config so it's not baked into the graph. In that case we'd always find an adjustment factor when building the graph, but only add the adjustment constant when preparing the responses if geoidElevation=true in the router config. Either way I think we want to include exactly the same new fields in the response: the constant that must be added to the ellipsoid to get the geoid, and a boolean for whether the elevations are relative to the geoid rather than the ellipsoid.

ellipsoidToGeoidDifference: 40,
geoidElevation: true (would be false if elevation figures are relative to the ellipsoid)

[barbeau]

Also FYI, I'd definitely suggest writing a number of unit tests to make sure EarthGravitationalModel is giving the proper output. I did a port of this class to Android, and while it passes the single unit test from GeoTools it didn't yield correct results when trying other values. Something definitely could have gone wrong with my port, I haven't had time to investigate yet. More info, including new unit tests that are failing, are at barbeau/earth-gravitational-model#1.

[demory]

We discussed this further offline today and think we have an agreement to proceed using the approach @abyrd suggests two posts above, with the ellipsoid/geoid difference computed for each graph and returned with a response, and a boolean configuration option to apply this to the returned elevation values. This flag will also be included with the response allowing the client to construct both elevation values. Also, per @barbeau's suggestion I will also confirm that EarthGravitationModel is producing the expected results.

This will address the immediate concern raised by TriMet, though we should note there are other options we could look at longer-term -- in particular, supporting Mapzen's elevation service (https://mapzen.com/documentation/terrain-tiles/formats/) as an alternative to using NED data directly.

[demory]

Hi all -- I've completed a basic implementation of the above which I've pushed to the branch geoidElevation.

Note that in the course of writing/testing this, I discovered that there are several different geoids that have been associated with WGS84. As best I can tell, the egm180.nor data (which the GeoTools EarthGravitationalModel is based on) refers to the EGM84 geoid, rather than the newer EGM96 and EGM2008 models which some of the above utilities use. Here is a handy tool for calculating all three at once: http://geographiclib.sourceforge.net/cgi-bin/GeoidEval

I tested EarthGravitationalModel as used in this OTP branch for about a dozen places around the world and got exact matches with EGM84 values from that calculator. @barbeau -- the unit tests you referenced appear to be written for EGM96, which is perhaps why you saw some failures in your port of the Geotools class. Also, those tests appear to use (0-360) longitudes, whereas I am using (-180-180).

@fpurcell and @grant-humphries -- According to EGM84 the ellipsoid/geoid difference for downtown Portland is -19.2m, which is what the modified OTP code will use. @abyrd previously cited a figure of about -22m, which looks correct for the EGM96 geoid. A number of ~40m was cited at the beginning of this thread; I'm not sure where that came from, but I'm pretty confident at this point that the correct offset is in the 19-22m range. I am about to rebuild/redeploy our Portland test server using this branch -- will let you know when it's ready.

[grant-humphries]

@demory with respect to the 40m figure that I initially cited I think I just quickly did a google search and found an online tool to get values to compare to what OTP was outputting, so the comparison was quickly thrown together. I defer to the work you're doing which looks to be a lot more thorough and sound.

[abyrd]

Looks like we got this geoid-ellipsoid difference problem completely sorted out back in 2017. Closing.