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

score 1 · Accepted Answer · answered Jul 27 '22 at 15:22

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

Thank you very much for your answer. Now the length is nearly the same as that we get from ArcGIS pro. But still, there is a difference of 0.601665 meters. Can we somehow go close to the length that we get from ArcGIS Pro using some other approach? — Rajesh Patnaik, Jul 28 '22 at 07:39
0.60 meters is less than a 0.1% difference and likely as good as you can get. — Shawn, Jul 28 '22 at 14:03

Calculate length using GeoPandas

1 Answers1