<jupyter_start><jupyter_text>Compare Geographic Distance Calculation Approaches<jupyter_code>%load_ext lab_black
%load_ext autoreload
%autoreload 2
import math

import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import Point<jupyter_output><empty_output><jupyter_text>[Table of Contents](table-of-contents)0. [About](about)1. [User Inputs](user-inputs)2. [Generate Data](generate-data)3. [Get Haversine Distance](get-haversine-distance)4. [Get Distance using `distance` method of `geopandas` library](get-distance-using-`distance`-method-of-`geopandas`-library)5. [Get Difference in Distances](get-difference-in-distances)  0. [About](about) This is a comparison of two approaches to calculate the distance (in kilometers) between two locations based on their latitude and longitude.  1. [User Inputs](user-inputs)<jupyter_code>l = 1_000_000
crs = 4326
epsg = 6345

# Bounds for city of Chicago from visual inspection of city on Google Maps
lat_extrema = [-87.7, -87.518]
lon_extrema = [41.64, 42.03]<jupyter_output><empty_output><jupyter_text>2. [Generate Data](generate-data) We'll generate one million random latitudes and longitudes corresponding to locations within the bounds of the city of Chicago. These bounds were picked at some points located outside the [city on Google Maps](https://www.google.com/maps/place/Chicago,+IL,+USA/@41.8339037,-87.8720459,11z). We'll generate the lat-long coordinates within each extreme pair of latituds and similarly for longitudes.Since we need to compare these locations to other locations, we'll generate two such datasets each with one million random locations' latitude-longitude coordinates.<jupyter_code>%%time
df1 = pd.DataFrame(
    zip(*[np.random.uniform(lon_extrema[0], lon_extrema[1], l), np.random.uniform(lat_extrema[0], lat_extrema[1], l)]),
    columns=["lat", "lon"],
)
df2 = pd.DataFrame(
    zip(*[np.random.uniform(lon_extrema[0], lon_extrema[1], l), np.random.uniform(lat_extrema[0], lat_extrema[1], l)]),
    columns=["lat", "lon"],
)
df = df1.assign(lat_end=df2["lat"]).assign(lon_end=df2["lon"])
display(df.head(2))<jupyter_output><empty_output><jupyter_text>3. [Get Haversine Distance](get-haversine-distance) Calculate the distance [using](https://stackoverflow.com/a/53490513/4057186) the [haversine formula](https://en.wikipedia.org/wiki/Haversine_formula) as the [estimation of the shortest point-to-point distance on a sphere](https://en.wikipedia.org/wiki/Great-circle_distance)<jupyter_code>def haversine(s_lat, s_lng, e_lat, e_lng):

    # approximate radius of the Earth in km
    R = 6373.0

    s_lat = s_lat * np.pi / 180.0
    s_lng = np.deg2rad(s_lng)
    e_lat = np.deg2rad(e_lat)
    e_lng = np.deg2rad(e_lng)

    d = (
        np.sin((e_lat - s_lat) / 2) ** 2
        + np.cos(s_lat) * np.cos(e_lat) * np.sin((e_lng - s_lng) / 2) ** 2
    )

    return 2 * R * np.arcsin(np.sqrt(d))
%%time
df["haversine"] = pd.Series(
    np.vectorize(haversine)(
        *[
            df[c]
            for c in ["lat", "lon", "lat_end", "lon_end"]
        ],
    ), index=df.index
)
display(df.head())<jupyter_output><empty_output><jupyter_text>4. [Get Distance using `distance` method of `geopandas` library](get-distance-using-`distance`-method-of-`geopandas`-library) Calculate the distance using the Python [`geopandas` library's `distance()` method](https://geopandas.readthedocs.io/en/latest/docs/reference/api/geopandas.GeoSeries.distance.html). This can be called on a `GeoSeries` (a column of [`Shapely`](https://pypi.org/project/Shapely/) data objects such as `Point`, `geometry`, or others - see [here](https://autogis-site.readthedocs.io/en/latest/notebooks/L1/geometric-objects.html)) with the distance to be computed between that column and another `GeoSeries`. First, we'll create two `GeoDataFrame`s - one with each lat-long coordinate from a single dataset generated earlier<jupyter_code>%%time
gdf_start = gpd.GeoDataFrame(
    df1[["lat", "lon"]],
    geometry=gpd.points_from_xy(df1["lon"], df1["lat"]),
    crs=crs,
)
gdf_end = gpd.GeoDataFrame(
    df2[["lat", "lon"]],
    geometry=gpd.points_from_xy(df2["lon"], df2["lat"]),
    crs=crs,
)
display(gdf_start.head(2))
display(gdf_end.head(2))<jupyter_output><empty_output><jupyter_text>We'll now calculate the distance between each `GeoDataFrame`'s `GeoSeries` (the `geometry` column)<jupyter_code>%%time
df["geopandas"] = (
    gdf_end["geometry"]
    .to_crs(epsg=epsg)
    .distance(gdf_start["geometry"].to_crs(epsg=epsg))
) / 1_000<jupyter_output>CPU times: user 4.58 s, sys: 52.7 ms, total: 4.63 s
Wall time: 4.64 s<jupyter_text>**Notes**1. In order to calculate this distance, we have to set an initial Coordinate Reference System (CRS) when defining the `GeoDataFrame`. The CRS indicates how the specified point coordinates are related to places on the Earth. In the above code, this is done by specifying the `crs` keyword argument to the `GeoDataFrame` class. Then, we need to project these coordinates into a local coordinate system ([UTM](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system) - a rectilinear projection coordinate system) where distance can be calculated with pythagorean theorem with good precision. Keeping in mind that we want to get the distance in meters, the CRS used for the city of Chicago was taken from [this ArcGIS 10.4.0 Projected Coordinate System lookup table](https://desktop.arcgis.com/en/arcmap/10.4/map/projections/pdf/projected_coordinate_systems.pdfpage=60) (see page 60 for the range of latitudes and longitudes that encompass the city of Chicago) and can be visualized [here](https://epsg.io/6345). This CRS conversion was made by calling the `to_crs()` method in the code above.  5. [Get Difference in Distances](get-difference-in-distances) We'll now calculate the absolute and percent difference between the two approaches relative to the `haversine` distance<jupyter_code>df["diff_km"] = df["geopandas"] - df["haversine"]
df["diff_pct"] = (df["diff_km"] / df["haversine"]) * 100
display(df.head())
display(df[["diff_km", "diff_pct"]].describe())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12940026075619296, 'fraction_numerical': 0.025097783572359842, 'mean_word_length': 5.236788617886178, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 10, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '252836', 'license': 'MIT', 'n_tokens_mistral': 2195, 'n_tokens_neox': 2059, 'n_words': 560, 'path': '1_3_compare_distance_calculation_approaches.ipynb', 'repo_name': 'edesz/chicago-bikeshare', 'size': '6136'}