Calculate the distance between two coordinates (WGS84) in ETRS89

Question

I have an Excel-file with coordinates (WGS84). From these coordinates I create a line (shp) and transform it (ETRS89). After the transformation I calculate the distance from point to point. I do this using a Python-script and ArcGIS.

However, I would like to calculate the distance between the points before I run my script and without an ArcGIS license. I found several answers to that question (i.e. John Cook,Heversine I,Heversine II and Heversine III). But they are all not accurate enough. I would like to get exact the same value as I get using "Calculate Geometry" in ArcGIS in ETRS89) without using ArcGIS. There is only a difference of a few centimeters, but in sum it is a few kilometers. Is there any option to do that? Any special libaries, which allow to transform coordinates and calculate the distance?

Answer 1

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]]

Answer 2

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,

Answer 3

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

Answer 4

from pyproj import Geod
g = Geod(ellps='WGS84')
az12,az21,dist = g.inv(-12.0928828908918,-77.06514436070439, -12.09163432452833,-77.06410713903763)
print(dist)

119.92574021957807