GeoPandas combines the capabilities of pandas and Shapely, making it easy to perform spatial operations. One of its most useful features is the ability to calculate distances between geographic points. In this tutorial, we’ll explore how to use GeoPandas to calculate distances using a simple dataset of city locations. We’ll cover how to compute distances between points, create distance matrices, and find the closest pairs of cities in the dataset. 

Setting Up Our Environment

First, let’s set up our environment and create our sample dataset:

# Import necessary libraries
import geopandas as gpd
from shapely.geometry import Point

# Create sample data
data = {'city': ['New York', 'Los Angeles', 'Chicago', 'Houston'],
        'geometry': [Point(-74.006, 40.7128), Point(-118.2437, 34.0522), 
                     Point(-87.6298, 41.8781), Point(-95.3698, 29.7604)]}
gdf = gpd.GeoDataFrame(data, crs="EPSG:4326")

print("Initial GeoDataFrame:")
print(gdf)

This code snippet creates a GeoDataFrame containing four major U.S. cities and their corresponding latitude and longitude coordinates. We’re using the EPSG:4326 Coordinate Reference System (CRS), which is the standard for geographic coordinate systems using latitude and longitude.

Expected output:

Initial GeoDataFrame:

          city                geometry
0     New York  POINT (-74.006 40.713)
1  Los Angeles  POINT (-118.24 34.052)
2      Chicago   POINT (-87.63 41.878)
3      Houston    POINT (-95.37 29.76)

In the following sections, we will explore few distance calculation scenarios based on our sample GeoDataFrame. 

1. Calculating distance from one point to all others

Often, you might want to calculate the distance from one specific point to all other points in your dataset. Here’s how you can do that:

# Project to a suitable coordinate system (e.g., USA Contiguous Albers Equal Area Conic)
gdf_projected = gdf.to_crs('ESRI:102003')

# Calculate distance between New York and all other cities
ny = gdf_projected[gdf_projected['city'] == 'New York'].geometry.iloc[0]
gdf_projected['distance_to_ny_km'] = (gdf_projected.geometry.distance(ny) / 1000).round(2)

print("\nGeoDataFrame with distances to New York (in km):")
print(gdf_projected[['city', 'distance_to_ny_km']])

In this example, we first project our data to a coordinate system suitable for distance calculations across the contiguous United States. Specifically, we use the ESRI:102003 projection, which is the USA Contiguous Albers Equal Area Conic projection. This projection is ideal for our purposes because:

1. It’s designed specifically for the contiguous United States, making it appropriate for our dataset of U.S. cities.
2. It’s an equal-area projection, which means it preserves area measurements across the map.
3. It uses meters as its unit of measurement, which allows for accurate distance calculations.

After projecting the data, we select the geometry of New York City and use the distance method to calculate the distance from New York to all cities in the GeoDataFrame. The results are stored in a new column called ‘distance_to_ny_km’. 

Expected output:

GeoDataFrame with distances to New York (in km):
          city  distance_to_ny_km
0     New York               0.00
1  Los Angeles            3912.96
2      Chicago            1138.57
3      Houston            2275.75

This approach ensures that our distance calculations are as accurate as possible for the geographic area we’re working with.

2. Creating a distance matrix

Sometimes, you might need to calculate distances between all pairs of points in your dataset. This is where a distance matrix comes in handy:

# Calculate distance matrix with city names as labels
distance_matrix = gdf_projected.geometry.apply(lambda g: gdf_projected.geometry.distance(g) / 1000).round(2)
distance_matrix.index = gdf_projected['city']
distance_matrix.columns = gdf_projected['city']
print("\nDistance Matrix (in km):")
print(distance_matrix)

This code creates a matrix where each cell represents the distance between two cities in kilometers. The rows and columns are labeled with city names for easy interpretation. The diagonal of this matrix will always be zero, as it represents the distance from a city to itself.

Expected output:

Distance Matrix (in km):
city         New York  Los Angeles  Chicago  Houston
city                                                
New York         0.00      3912.96  1138.57  2275.75
Los Angeles   3912.96         0.00  2790.06  2200.12
Chicago       1138.57      2790.06     0.00  1520.39
Houston       2275.75      2200.12  1520.39     0.00

3. Finding the closest pair of points

Once you have a distance matrix, you can use it to find the pair of points that are closest to each other:

# Find the two closest cities
min_distance = distance_matrix.where(distance_matrix > 0).min().min()
closest_pair = distance_matrix[distance_matrix == min_distance].stack().index.tolist()[0]
print(f"\nThe two closest cities are {closest_pair[0]} and {closest_pair[1]}")
print(f"They are {min_distance} km apart.")

This code first finds the minimum non-zero distance in the matrix (we use where(distance_matrix > 0) to ignore the zero distances along the diagonal). Then, it identifies the pair of cities that are separated by this minimum distance. By using city names as labels in our distance matrix, the output of this code is now more intuitive and directly references the city names.

Expected output:

The two closest cities are New York and Chicago
They are 1138.57 km apart.

Conclusion

In this tutorial, we’ve explored how to use GeoPandas to perform various distance calculations on geospatial data, including calculating distances from one point to all others, creating a distance matrix, and finding the closest pair of points. It’s worth noting that while the ESRI:102003 projection is excellent for the contiguous United States, you might need to choose a different projection for datasets covering other parts of the world or for global datasets.