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.

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

|improve this answer|||||

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.