Calculate length using GeoPandas

Question

I have a shape file "Line1.shp". It has only one feature. When I calculate the length in ArcGIS Pro using "Calculate Geometry" and selecting units as "Meters" and coordinate system as "GCS_WGS_1984" it gives 1919.182971

I am trying to calculate the length of line features (in meter unit) using GeoPandas. I use the following code:

import geopandas as gpd
fn = "D:/Line1.shp"
df = gpd.read_file(fn)
df = df.to_crs('EPSG:3857')
ps_length = df['geometry'].length
print(ps_length)

It returns:

**>>> print(ps_length)

0 2591.493383

dtype: float64 **

Using GeoPandas I am getting 2591.493383 as length. I think the above code will return the length in meters units as the unit of measurement of 'EPSG:3857' is in meters.

Why am I getting a different length using these codes? Am I doing anything wrong in choosing the crs?

I have uploaded the shape files in the link. https://drive.google.com/file/d/1TkibJLQJatn2bzNUSM2nRP7TFdYkw3bq/view?usp=sharing

Answer 1

The 3857 CRS can have inaccurate length measurements away from the equator. See the answer here: How can EPSG:3857 be in meters?

ArcGIS is likely converting to the appropriate UTM CRS instead. Here is some code to do that in python. I've converted your line to a WKT string for simplicity.

import shapely
import geopandas as gpd
from pyproj import CRS
import utm


line_wkt = 'LINESTRING (-79.99193914061688 42.17279348571836,' \
           '-79.99182484706496 42.172693810334636,'\
           '-79.99174384148422 42.172616260812156,'\
           '-79.98826035896703 42.169181672722615,'\
           '-79.98780764938634 42.16871996165105,'\
           '-79.98759224056774 42.16849920349132,'\
           '-79.98734240740913 42.16834931134773,'\
           '-79.98708048113015 42.1682301120582,'\
           '-79.9869009456341 42.16810948504351,'\
           '-79.98673382755022 42.16789181879166,'\
           '-79.98652224016483 42.167552422453554,'\
           '-79.98621147967039 42.16718129057849,'\
           '-79.9823844319812 42.16335683544205,'\
           '-79.9813640545146 42.16233704579152,'\
           '-79.97830098894026 42.15927548394874,'\
           '-79.97804026370885 42.158999913955604)'


def utm_crs_from_latlon(lat, lon):
    crs_params = dict(
        proj = 'utm',
        zone = utm.latlon_to_zone_number(lat, lon),
        south = lat < 0
        )
    return CRS.from_dict(crs_params)

line_shape = shapely.wkt.loads(line_wkt)

utm_crs = utm_crs_from_latlon(lat = line_shape.centroid.y,
                              lon = line_shape.centroid.x)

line_geo   = gpd.GeoSeries([line_shape], crs='EPSG:4326')

line_geo.to_crs('EPSG:3857').length
> 2591.493383
line_geo.to_crs(utm_crs).length
> 1918.581306