1

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?

  • 1
    github.com/jswhit/pyproj proj4 should work for you – Alex Tereshenkov Jan 26 '16 at 10:38
  • 2
    could you specify the distance that you are interested in. Basically, the length of a straight line joining two points is not the same as the shortest arc between those points. It is unclear to me what you want to measure. Also you should have the same length in ETRS 89 and in WGS84 – radouxju Jan 26 '16 at 11:14
  • @radouxju I thought I had specified that in my question: I would like to optain the length of the line after I transformed it to ETRS89 using the calculate geometry method in ArcGIS. For this optian, as far as I know, I cannot specify if its planar or geodesic, or can I? – Alex Jan 26 '16 at 11:34
  • 1
    All of the Haversine formulas use a sphere. To match that in ArcGIS, you'd have to have the data's CRS use the same sphere model. With lat/lon data, ArcGIS is using a geodesic calculation (roughly Vincenty). You're not going to be able to match it even by adjusting the sphere radius in a Haversine formula. – mkennedy Jan 26 '16 at 17:12
  • 3
    @Alex Could you please provide more information about your coordinate system. You mention ETRS89, which is a geographic CRS, but as far as I know the calculate geometry tool should be disabled in ArcGIS if your line is in Lat/long. Therefore my guess is that you are not working in ETRS89 but you use a projected coordinate system that is based on ETRS89. In this case, you should not expect your line to have the same length as the geodesic distance that you approximate with haversine. If you want the same value as in ArcGIS calculate geometry, use the euclidian distance in the projected space. – radouxju Jan 26 '16 at 19:14
3

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

  • Yes you are correct. I was not aware that ETRS89 is also an geographic coordinate system, as I always use WGS84 as a geographic cs and ETRS89/UTM as a projected coordinate sytem. Now much is clearer to me and I actually would like to calculate the "!shape.geodesicLength@meters!"-value without using ArcGIS. Is there a library to do that? The gdaltransform link looks really interesting. However, I don't really get how to imply that to my python script. Could you give some more details on this? – Alex Jan 27 '16 at 8:57
  • in this case, doesn't this post help gis.stackexchange.com/questions/163785/… ? – radouxju Jan 27 '16 at 9:33
  • and if you want the most accurate, google Vincenty – radouxju Jan 27 '16 at 9:34
  • No, the post you mention describes the heversine function, which is not accurate enough. The Vincenty stuff looks nice, but I tried for instance the distance calculation using vincenty from geopy, which gives values not accurate enough, too. Well, accuracy in relation to the calculated !shape.geodesicLength@meters! in ArcGIS. – Alex Jan 27 '16 at 9:42
1

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,

  • I get closer that way, but there is still a difference. – Alex Jan 26 '16 at 11:42
  • Well, there is still a difference of about 10 cm for a distance of 50 m. As the total distance is about 700 km, the calculations differ of about 1.6 km. Using QGIS I have get a difference of 0.1 cm for a distance of 50 m and a total difference of 35 m. I could live with the QGIS-value, but a km is definetly too much. – Alex Jan 26 '16 at 12:55
  • @Alex, from doing some checking on the internet, it seems that different formulas have pros and cons under different conditions. e.g. over very big distances, or very small distances. It seems you are using lots of small distances, so the Haversine formula should be ok. (See movable-type.co.uk/scripts/latlong.html) – Reversed Engineer Jan 27 '16 at 10:43
  • Another thought - I'm not really familiar with Python, but is there any way to control numerical precision? E.g. many languages there is single and double precision numerical types (e.g. in Pascal and C# etc) there is the single and double numerical types. Does Python have something similar? You'd be better off using the highest precision floating point possible – Reversed Engineer Jan 27 '16 at 10:46
  • Are all the errors out by the same percentage (up or down), or are they more or less the same? – Reversed Engineer Jan 27 '16 at 10:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.