Calculate the distance between two coordinates (WGS84) in ETRS89

The error is no doubt because the earth radius in your formula is not exactly the same as that used by ArcGIS. In fact, seeing as the earth is not a perfect sphere, the radius is different at the equator as it is at the poles. Probably ArcGIS corrects for that.

However: in the Haversine Python script, it has:

Base = 6371 * c

If you calculate the apparent radius of the earth where you are (or your map is), and adjust the radius in the code, it will be much more accurate. I would just determine the percentage difference between what you expect and what you get, and adjust the 6371 kilometres in the code above accordingly.

Getting exactly the same distance (what? to within one millimeter? one micrometer?) is just a dream. Don't waste time trying, unless you have the exact formula used by ArcGIS. I'm sure it's compensating for the non-spherical shape of the earth, and maybe even for differences in altitude as well,

WGS 84 and ETRS 89 are two geographic coordinate systems (Lat/long). With those coordinate system, you will measure distances on the surface of the ellipsoid. WGS84 and ETRS 89 use almost identical spheroid (see below), so in most cases you will not see any difference between the 2.

You are projecting your data in Universal Transverse Mercator zone 35 (based on ETRS 89 datum). UTM projection is conformal, so it preserves angles and approximates shape but distorts distance and area. This means that the length of your segment between two points projected in UTM will not be exactly the same as the geodetic distance between those points.

In practice, if you want to get the same length as the result of a ArcGIS "calculate geometry" field calculation, you should project your points (e.g. gdaltransform -a_srs EPSG:4326 -t_srs EPSG:25832 sourcefile outputfile), then you compute the euclidian distance between your points (sqrt((x_a-x_b)²+(y_a-y_b)²))

Finally, with recent ArcGIS versions, you can also compute the geodetic length by using the following command in the field calculator :

!shape.geodesicLength@meters!

GEOGCS["ETRS89",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

def wgs84_dist(lat1, lon1, lat2, lon2):
    """
    Calculates distance based on WGS84 Geode
     lat1, lon1 = origin
     lat2, lon2 = destination
    """

    if lat1 > 1000:
        print('[fs-utils] - ERROR: Wrong input. Needs to be decimal degrees ')

    dist = geod.Inverse(lat1, lon1, lat2, lon2)['s12']

    return dist