5

I need to calculate the distance between points of two different datasets, once with locations of towns, the other with locations of dams. Both are located in different countries in Africa.

The more I read the more confused I keep getting. I'm using the Python ecosystem (Geopandas, Shapely, Fiona).

A) The first dataset is a .shp file with locations of towns as Points. This one is nice and provides the crs as epsg=4326 (which I understand is the code for WGS 84):

import geopandas as gpd
towns = gpd.read_file('towns.shp')) #Uses fiona to load
print(towns.crs)
print(towns.geometry.head(3))
#{'init': 'epsg:4326'}
#0    POINT (8.877318824300001 9.93427297769)
#1    POINT (9.163896418389999 9.47532028526)

B) The second one has locations of dams. It came in an excel file with latitude and longitude with decimal places. It didn't give any information on the crs.

#Omitting the excel loading of `dams`
print(dams.geometry.head(2))
print(dams.crs) # empty
#487    POINT (8.97333333333 9.76472222222)
#488    POINT (4.55305555556 8.44277777778)

What's the right way to measure the distance between these two points? (Let's assume for now they are both in WGS 84)

I can calculate the distance between all towns and the first element of the dams:

# Distance function uses Shapely
print(towns.geometry.distance(dams.geometry.iloc[1]))
#0       0.194849
#1       0.346508
#2       1.046174

Is this just calculating the Euclidian distance and hence very inaccurate?

What can I do for this workflow?

Should I transform the crs of the points to something that works for all of Africa and then take the Geopandas/Shapely distance function or would it be easier to keep the lat/lon (or WGS 84) and use a Haversine formula (or similar)?

To me this would break a bit the benefit of using Geopandas.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
cd98
  • 151
  • 1
  • 6
  • 1
    Welcome to GIS SE. As a new user, please take the Tour. It is general policy here to have exactly one question per Question. Please edit your question to focus on your more important topic (Given the coarse precision -- 3 places is ~111m -- it probably wouldn't hurt to assume WGS84; only the provider can tell you for sure). – Vince Nov 25 '15 at 16:48
  • @Vince: thanks! I think that answers my first question. I've edited to include just the second (and main) one. – cd98 Nov 25 '15 at 17:14

1 Answers1

2

1. Original trigonometric calculation

I have been using a direct trigonometric calculation function from before I got on python, had got this from a SE forum answer. It works directly with lat-long values.

from math import sin, cos, sqrt, atan2, radians
def lat_long_dist(lat1,lon1,lat2,lon2):
    # function for calculating ground distance between two lat-long locations
    R = 6373.0 # approximate radius of earth in km. 

    lat1 = radians( float(lat1) )
    lon1 = radians( float(lon1) )
    lat2 = radians( float(lat2) )
    lon2 = radians( float(lon2) )

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))

    distance = float(format( R * c , '.2f' )) #rounding. From https://stackoverflow.com/a/28142318/4355695
    return distance

# test:
lat_long_dist(17.82149,78.331,17.34932,78.54818)

>>> 57.35

2. Projection change

When using QGIS, I found that I had to convert all my shapefiles to EPSG:3857 (pseudo mercator) or another projection that converts the lat-longs to another set of numbers (my guess is meters distance from equator), only then I could do constant distance buffer, grid, etc in meters. See this post: create buffer in meters. There's a caveat that EPSG:3857 gets inaccurate nearer the poles. My locations were nearer to equator so it didn't matter in my case.

So even in your towns.geometry.distance(point) case I'll reckon you have to convert your data to a CRS. Since your data is in Africa, EPSG:3857 should be ok.

Nikhil VJ
  • 1,026
  • 12
  • 31