3

I am having problems getting lengths from multilinestrings after transforming data from EPSG:27700 - OSGB 1936 / British National Grid to EPSG:4326/ WGS 84.

I'm doing this conversion to import a roads shapefile into a PostGIS database that uses WGS84 as default. However, when I try to calulate lengths, I'm having weird readings.

Firstly I thought about some incorrect transformation, but after doing some additional tests in an empty database importing the layer without transforming it (OSGB 1936 reference), I´m still puzzled:

The following query returns contradictory results

SELECT  
    ST_Length(geom) AS length_OSGB
    ,ST_Length_Spheroid(ST_Transform(geom,4326),'SPHEROID["WGS 84",6378137,298.257223563]') As length_WGS84    
FROM osgb36_data
WHERE id = 108   

this yields 338 meters for the OSGB36 length and 482 meters for the WGS84 length. The right one is 338 meters.

I've done the importation to the database using QGIS and DBmanager.

Also, I had to guess the initial projectios, as the .prj was missing, but projecting the result over google maps results in a perfect match.

Any hint is welcome

EDIT:

This is the WKT for the element in the example, as suggested:

    SRID=27700;MULTILINESTRING((
 423216.279 574665.249 0,
 423206.315 574708.077 44.649,
 423158.458 574896.911 242.525,
 423132.406 574993.061 344))

Well, that's weird, seems like they were using the z coordinate to register the accumulated distance for each point in the linestring... That could explain the issue with the wrong distances

4
  • 1
    Paste the WKT of the object itself (ST_AsEWKT(geom)) so others can verify. Commented Jul 7, 2015 at 15:20
  • How do you know that 338m is the right answer?
    – mlinth
    Commented Jul 7, 2015 at 19:30
  • One of the fields of the feature is the total length of the segment. I'll post the wkt tomorrow . Thanks for the interest Commented Jul 7, 2015 at 19:32
  • ewkt posted...! Commented Jul 8, 2015 at 5:38

2 Answers 2

4

The length functions work differently with 3D linestring geometries:

  • ST_Length - returns 2D distances for geometry types, and oddly 3D distances for geography types (but not in this question)
  • ST_Length_Spheroid - returns 3D distances for geometry types

Your example is in 3D, so it will calculate the 3D length with ST_Length_Spheroid and the 2D length with ST_Length on a geometry.

However, if you always want 2D lengths, you can force the geometry to a 2D linestring using either ST_Force2D or ST_Force_2D function (the name change in PostGIS 2.1).

SELECT ST_Length(geom) AS length_27700,
  ST_Length_Spheroid(ST_Transform(ST_Force2d(geom), 4326), spheroid) AS length_spheroid_2d,
  ST_Length_Spheroid(ST_Transform(geom, 4326), spheroid) AS length_spheroid_3d
FROM (
  SELECT 'SRID=27700;MULTILINESTRING((
      423216.279 574665.249 0,
      423206.315 574708.077 44.649,
      423158.458 574896.911 242.525,
      423132.406 574993.061 344))'::geometry AS geom,
    'SPHEROID["WGS 84",6378137,298.257223563]'::spheroid
) AS f;

  length_27700   | length_spheroid_2d | length_spheroid_3d
-----------------+--------------------+--------------------
 338.39264086563 |   338.515774984236 |   482.622196178051
3
  • so why is different ST_Length(geom) and ST_Length_Spheroid(ST_Transform(geom, 4326), spheroid)? Aren't both of them 3d coords? Commented Jul 8, 2015 at 9:56
  • ST_Length is documented to calculate a 2D distance for geometries. ST_Length_Spheroid is documented to behave differently.
    – Mike T
    Commented Jul 8, 2015 at 9:59
  • Thanks. I would add that clarification in the answer, that would help a lot. I'm marking this as the answer, although both yours and @nickves are good answers Commented Jul 8, 2015 at 10:03
1

Postgis is correct.

Your line is in 3d space, and in fact has a length of ~482m even in the British National Grid.

338 meters is also correct but it corespondents in the projection of the 3d line in the 2d space.

 with a as ( select st_Geomfromewkt('SRID=27700;MULTILINESTRING((
  423216.279 574665.249 0,
  423206.315 574708.077 44.649,
  423158.458 574896.911 242.525,
  423132.406 574993.061 344))') geom) 
   select ST_Length(geom) "2D Length", ST_3DLength(geom) "3D_Length" from a;


     2d_length    |    3d_length
 -----------------+-----------------
  338.39264086563 | 482.54086040695 (1 row)

Manual Links:
ST_Length
ST_3DLength

Edit:
Related:
http://revenant.ca/www/postgis/workshop/measurement.html

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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