Commit 9f4dbd8d authored by Fabian Kovac's avatar Fabian Kovac
Browse files

[i] updated distance function using a wgs84 ellipsoid

parent 4984ea71
......@@ -67,32 +67,42 @@ def _log(msg: str) -> None:
def get_distance(lat_a: pd.Series, lon_a: pd.Series, lat_b: pd.Series, lon_b: pd.Series) -> np.array:
"""Calculcates distance between two coordinates in km using the haversine function
"""Calculcates distance between two coordinates in km
using a rotation-ellipsoid in cartesian coordinates out of polar coordiantes
Parameters:
lat_a (pd.Series): Latitudes of point A
lon_a (pd.Series): Longitudes of point A
lat_b (pd.Series): Latitudes of point B
lon_b (pd.Series): Longitudes of point B
Returns:
km (np.array): Vector with distances in km (can directly be assigned a pandas column)
"""
# convert latitudes and longitudes to radians
lat_a, lon_a = np.radians(lat_a), np.radians(lon_a)
lat_b, lon_b = np.radians(lat_b), np.radians(lon_b)
# calculate differences of latitudes and longitudes
diff_lat = lat_a - lat_b
diff_lon = lon_a - lon_b
# calculate distances in km between point a and point b
a = np.sin(diff_lat/2)**2 + np.cos(lat_a) * np.cos(lat_b) * np.sin(diff_lon/2)**2
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
km = 6371 * c
return km
Returns:
length (np.array): Vector with distances in km (can directly be assigned a pandas column)
"""
# constants (euqator radius and pole radius in km)
# r_equator is the 'semi-major axis' and r_pole the 'semi-minor axis' on a WGS84 ellipsoid
r_equator = 6378.137
r_pole = 6356.7523142
# calculate rotation-ellipsoid in cartesian coordinates out of polar coordiantes
za = np.sin(np.radians(lat_a)) * r_pole
ra = np.cos(np.radians(lat_a)) * r_equator
xa = np.sin(np.radians(lon_a)) * ra
ya = np.cos(np.radians(lon_a)) * ra
zb = np.sin(np.radians(lat_b)) * r_pole
rb = np.cos(np.radians(lat_b)) * r_equator
xb = np.sin(np.radians(lon_b)) * rb
yb = np.cos(np.radians(lon_b)) * rb
# calculate distances between point a and point b
dx = xa - xb
dy = ya - yb
dz = za - zb
length = np.sqrt(np.square(dx) + np.square(dy) + np.square(dz))
return length
def prep(file_config: pathlib.Path, file_trans: pathlib.Path) -> None:
......@@ -158,9 +168,9 @@ def prep(file_config: pathlib.Path, file_trans: pathlib.Path) -> None:
_log('Removed duplicated links which are not in use')
# calculate LENGTH in km between links using the haversine function
# calculate LENGTH in km between links
df_config['LENGTH'] = get_distance(df_config['LATITUDE_A'], df_config['LONGITUDE_A'], df_config['LATITUDE_B'], df_config['LONGITUDE_B'])
_log('Calculate distance between sites using the haversine function')
_log('Calculated distances between sites using a WGS84 ellipsoid')
# convert FREQUENCY to float
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment