1

I want to know linestrings (epsg:4326) length in meters instead of degrees. I know that I can use proj4 for that.

Here is my code, if you interested.

project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:3395')
)

project2 = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:3395'),
    pyproj.Proj(init='EPSG:4326')
)


def distance_in_meters(linestring):
    return transform(project, linestring).length


def interpolate(distance, linestring):
    current_ratio = 0
    points = []
    step = distance / distance_in_meters(linestring)

    while current_ratio < 1:
        point = linestring.interpolate(current_ratio, normalized=True)
        current_ratio += step
        points.append(point)
    return points


distance_between_points = 249.67
total_distance_in_meters = distance_between_points * 5

class TestInterpolate(TestCase):
    def test_interpolate(self):
        with open("linestring.json", 'r') as f:
            s = shape(json.load(f))
            points = interpolate(10, s)
            l = LineString(points)
            self.assertEqual(distance_in_meters(l), total_distance_in_meters)

And here is linestring that I use for tests.

{
  "type": "LineString",
  "coordinates": [
    [
      37.6232398,
      55.7615482
    ],
    [
      37.6262633,
      55.7606701
    ],
    [
      37.6232398,
      55.7615482
    ],
    [
      37.6262633,
      55.7606701
    ],
    [
      37.6232398,
      55.7615482
    ],
    [
      37.6262633,
      55.7606701
    ]
  ]
}

What I've done is I've taken coordinates of two points from google maps and measured distance between them also in google maps. Than I've duplicated it several times creating a path that swings back and forth 5 times.

And the result of the test is

1864.394777794176 != 1248.35

Which is pretty far off from the actual distance of linestring. I suppose that the reason is - I'm not using right CRS for meters projection.

0

Not an answer to your question but a reference method to compare your results with:

SELECT ST_Length(the_geog) As length_spheroid,  ST_Length(the_geog,false) As length_sphere
FROM (SELECT ST_GeographyFromText(
'SRID=4326;LINESTRING ( 37.6232398 55.7615482, 37.6262633 55.7606701, 37.6232398 55.7615482, 37.6262633 55.7606701, 37.6232398 55.7615482, 37.6262633 55.7606701)') As the_geog)
As foo;

Check the results:

length_spheroid 
1067.52749423872

length_sphere
1064.36864138641

I also converted the data into EPSG:32637 with OpenJUMP and measured the distance: 1067.1980.

Conclusion: The distance that you measured from Google maps is not a good reference. Are you sure that you measured that distance five, not six times? With five times it would make 1040.291 which is much closer to results with EPSG:32637 and with geography.

  • Is there a CRS with meters units that is covering the whole globe, not small part of it? – user1685095 Nov 23 '15 at 13:58
  • No, meters mean a 2d-rectangle of some sort, while the globe is a sphere. Meter-based projections always have a limited area of validity. – Uffe Kousgaard Nov 23 '15 at 14:43
  • @UffeKousgaard thanks. I have one more question though. I am not trying to calculate length of a geometry with geopy distance function (vincenty method). The error now is about 30 meters (if using postgis as reference. This seems strange to me. Can you help with that? Like maybe postgis and geopy using different spheroids or something? – user1685095 Nov 23 '15 at 15:56

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.