1

I know that a similar question has been asked multiple times here (a great answer can be found here: How can I get the geometry length in meters?), but in contrast I want to get the length of the reversed edge in meters.

I know that this is quite easy to achieve using ST_Length(ST_Reverse(geometry)), but look at the following queries (including their results as SQL comment; do not care about the id 1484 - it's just for getting a unique edge):

SELECT ST_Length(edge_geometry) FROM bolzano_street_edges WHERE edge_id = 1484;
--0.0007624185038570938 (length in SRID units... not meters)
SELECT ST_Length(ST_Reverse(edge_geometry)) FROM bolzano_street_edges WHERE edge_id = 1484;
--0.0007624185038570938 (length in SRID units... not meters)

This one works (meaning delivers the same length for the edge and it's reversed counterpart) and delivers the length in SRID units. Since my SRID is already EPSG:4326 getting the length in meters is easy using the geography type of PostGIS, right? Let's see...

SELECT ST_Length(edge_geometry::geography, true) FROM bolzano_street_edges WHERE edge_id = 1484;
--84.58170640981749 (length in meters)
SELECT ST_Length(ST_Reverse(edge_geometry)::geography, true) FROM bolzano_street_edges WHERE edge_id = 1484;
--84.58170640981747 (length in meters.... with 747 at the end; not 749)

Hmm... now that's something unexpected (at least to me).
The length of the edge is not the same as the length of its reverse. I know that the difference is very small, but for me that's a problem since I want to find edges with the same length by using their hash code. Even with this small difference this won't work.

Another approach that might work, is to use ST_LengthSperoid instead of the geometry type. According to the accepted solution from the question linked at the beginning of my question this should give the same result...

SELECT ST_LengthSpheroid(edge_geometry, 'SPHEROID["WGS 84",6378137,298.257223563]') FROM bolzano_street_edges WHERE edge_id = 1484;
-- 84.58170640981679 (length in meters... however, a bit different to the above)
SELECT ST_LengthSpheroid(ST_Reverse(edge_geometry), 'SPHEROID["WGS 84",6378137,298.257223563]') FROM bolzano_street_edges WHERE edge_id = 1484;
-- 84.58170640981679 (length in meters... same result for the reversed edge.. great)

I do have several questions regarding my obervations:

  • Is there a reason why the length of an edge is not the same as the length of its reversed counterpart when using the geography type? Is this a rounding error (or even a bug) in PostGIS or is this something expected?
  • Why is the result of ST_LengthSpheroid not exactly the same as the one of ST_Length(geography)?
  • Which result is the most accurate one??
6
  • Yes, this is a rounding error, a slight difference in the result depending on order of operands.
  • Not sure what ST_LengthSpheroid() differs, but it's possible there's a slightly different code path still, as they originally had completely different code paths.
  • You're talking about a difference of 0.00000000000002, there's no "most accurate" answer. (A difference of 0.00001 degrees is about 1m, so you're looking at a difference of about 2 nanometers: molecule scale.)

As to why things might differ at all, the code involves a converging solution so I'm not really surprised that order of arguments can result in a loop or two difference in when convergence is hit. Probably most systems that provide "consistent" results are just rounding away any differences long before a user sees them, maybe we should be doing more of that.

  • 1
    That the code uses a converging solution explains almost everything. Thanks for the information. I solved my problem by switching the datatype from double precision to real and using trunc after 10 digits when inserting. I do not need the accuracy (so saving some memory even improves performance a tiny bit) so this works well :-) – Nikolaus Krismer Aug 12 '16 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.