Calculate length of LineString in meters using Python

user_13272706 0
Viewed 349

I know this question has been asked before, but I just cannot figure it out. I'm just trying to find the length in meters from a LineString that I've made from two points ('EPSG:4326'). In R I get around 311 meters which is about the same as google maps straight line distance give using this code:

library(purrr)
library(sf)
library(dplyr)

xy <- list(c(59.94661244955431, 10.72052329576658), c(59.948009419166226, 10.725362000881619))

lat <- xy %>%
  map(1) %>%
  unlist()

lon <- xy %>%
  map(2) %>%
  unlist()

st_as_sf(data.frame(lon, lat), coords=c("lon", "lat"), crs="EPSG:4326") %>%
  summarise(do_union = TRUE) %>%
  st_cast("LINESTRING") %>%
  st_length()

311.0174 [m]

When I do it in Python I instead get 621 even if I reproject to 'EPSG:3857'. What am I doing wrong here? Python code below:

from pyproj import Transformer
from shapely.geometry import LineString
import shapely.ops as sp_ops


latlon = [[59.94661244955431, 10.72052329576658], [59.948009419166226, 10.725362000881619]]
xy = [(lon, lat) for lat, lon in latlon]
line = LineString(xy)
transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True)
line_transformed = sp_ops.transform(transformer.transform, line)
line_transformed.length

Out[147]: 621.7406578605734
shapely geospatial python
2 Answers

EPSG:3857 is useless for any kind of computation. It is nominally in meters but, because of the brutal distortion of the Mercator projection, those meters don't really mean anything. Taking a guess from your name and your sample coordinates your are in Norway, meaning that some UTM projection is probably a good bet.

user_41770090

Update: As mentioned I need to use the right projection which in this case is EPSG:32632 (UTM 32N):

from pyproj import Transformer
from shapely.geometry import LineString
import shapely.ops as sp_ops

latlon = [[59.94661244955431, 10.72052329576658], [59.948009419166226, 10.725362000881619]]
xy = [(lon, lat) for lat, lon in latlon]
line = LineString(xy)
transformer = Transformer.from_crs('EPSG:4326', 'EPSG:32632', always_xy=True)
line_transformed = sp_ops.transform(transformer.transform, line)
line_transformed.length

Out[2]: 311.9284270277389
user_132727060
Related