How to split epsg 4326 linestring in parts consisting of given length in meters in python?

Question

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.

Answer 1

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

• Read this PostGIS document http://postgis.net/docs/ST_Length.html
• Take the Geography example "Return length of WGS 84 geography line" and insert your linestring into the query

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.