1

I want to know the length in meters of a LineString with coordinates in EPSG4326.

Based on this post, this is code that works for me:

from shapely.geometry import LineString
from shapely.ops import transform
from functools import partial
import pyproj

roadsegment = LineString([(5.318945751388698, 50.91861496144046), (5.31929104222258, 50.91834578805272), (5.319299630296815, 50.91833937146759), (5.319523911140832, 50.918171650963316), (5.319554872486604, 50.91814850228405), (5.3195882431363435, 50.918123536051475), (5.319647895985286, 50.91807892920322), (5.319860389778624, 50.91792339522707), (5.320361197231897, 50.91757273001555), (5.320844582471854, 50.91724481403615), (5.3213449709422145, 50.916919032387966), (5.3217234583154465, 50.9166827732894), (5.322094602885187, 50.91646186473597), (5.32282743106912, 50.91603652941753)])

# Geometry transform function based on pyproj.transform
project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:32633'))

line2 = transform(project, roadsegment)
print(str(line2.length) + " meters")
#398.5297511214621 meters

It gives me 398.5297511214621 meters, which I have checked to be correct. However, it gives a FutureWarning about deprecated Proj initialization. In this post, I read that I could remove the ‘init=’ for the updated syntax. 

from shapely.geometry import LineString
from shapely.ops import transform
from functools import partial
import pyproj

roadsegment = LineString([(5.318945751388698, 50.91861496144046), (5.31929104222258, 50.91834578805272), (5.319299630296815, 50.91833937146759), (5.319523911140832, 50.918171650963316), (5.319554872486604, 50.91814850228405), (5.3195882431363435, 50.918123536051475), (5.319647895985286, 50.91807892920322), (5.319860389778624, 50.91792339522707), (5.320361197231897, 50.91757273001555), (5.320844582471854, 50.91724481403615), (5.3213449709422145, 50.916919032387966), (5.3217234583154465, 50.9166827732894), (5.322094602885187, 50.91646186473597), (5.32282743106912, 50.91603652941753)])

# Geometry transform function based on pyproj.transform
project = partial(
    pyproj.transform,
    pyproj.Proj('EPSG:4326'),
    pyproj.Proj('EPSG:32633'))

line2 = transform(project, roadsegment)
print(str(line2.length) + " meters")
# 636.8152519720558 meters

This indeed does not result in a warning anymore, but the outcome is now 636.82m. What I am doing wrong? I am new to working with geospatial data, so any help is welcome.

Marcelo Villa
  • 5,928
  • 2
  • 19
  • 38
tfdgrand
  • 13
  • 2

1 Answers1

2

The issue you are having is due to axis order changes in PROJ 6+: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6

One solution is to use always_xy: 

from shapely.geometry import LineString
from shapely.ops import transform
from functools import partial
import pyproj

roadsegment = LineString([(5.318945751388698, 50.91861496144046), (5.31929104222258, 50.91834578805272), (5.319299630296815, 50.91833937146759), (5.319523911140832, 50.918171650963316), (5.319554872486604, 50.91814850228405), (5.3195882431363435, 50.918123536051475), (5.319647895985286, 50.91807892920322), (5.319860389778624, 50.91792339522707), (5.320361197231897, 50.91757273001555), (5.320844582471854, 50.91724481403615), (5.3213449709422145, 50.916919032387966), (5.3217234583154465, 50.9166827732894), (5.322094602885187, 50.91646186473597), (5.32282743106912, 50.91603652941753)])

# Geometry transform function based on pyproj.transform
transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
line2 = transform(transformer.transform, roadsegment)
print(str(line2.length) + " meters")

The output is as expected: 398.5297511214502 meters

You can also get the geodesic length using pyproj: https://pyproj4.github.io/pyproj/stable/examples.html#geodesic-calculations

crs = pyproj.CRS("EPSG:4326")
geod = crs.get_geod()
length = geod.geometry_length(roadsegment)

Output is: 396.43737564550355

snowman2
  • 7,321
  • 12
  • 29
  • 54