Global static map

In this tutorial we will create a global map showing the connections between major airports around the world. Among other libraries we use in this tutorial, we will work with Cartopy to visualize our geospatial data. Cartopy is a Python library designed for cartography and geographic data processing. It’s built on top of Matplotlib and provides a powerful interface for creating maps and visualizing geospatial data. With Cartopy, you can easily manipulate and project geographic datasets onto different map projections, draw maps for different regions of the world, and overlay data points, lines, and polygons. Learn more about Cartopy here.

Data

In this tutorial we will work with the airports data from Openflights.org. This includes a rich database of airports around the world as well as information such as the geographical locations and global connections.

Let’s get started!

import pathlib
NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"

import pandas as pd
from shapely.geometry import Point
import geopandas as gpd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

airports = pd.read_csv(DATA_DIRECTORY / "airports.dat", delimiter=',',
                       names=['id', 'name', 'city', 'country', 'iata',
                              'icao', 'lat', 'long', 'altitude', 'timezone',
                              'dst', 'tz', 'type', 'source'])
airports.head()

	id	name	city	country	iata	icao	lat	long	altitude	timezone	dst	tz	type	source
0	1	Goroka Airport	Goroka	Papua New Guinea	GKA	AYGA	-6.081690	145.391998	5282	10	U	Pacific/Port_Moresby	airport	OurAirports
1	2	Madang Airport	Madang	Papua New Guinea	MAG	AYMD	-5.207080	145.789001	20	10	U	Pacific/Port_Moresby	airport	OurAirports
2	3	Mount Hagen Kagamuga Airport	Mount Hagen	Papua New Guinea	HGU	AYMH	-5.826790	144.296005	5388	10	U	Pacific/Port_Moresby	airport	OurAirports
3	4	Nadzab Airport	Nadzab	Papua New Guinea	LAE	AYNZ	-6.569803	146.725977	239	10	U	Pacific/Port_Moresby	airport	OurAirports
4	5	Port Moresby Jacksons International Airport	Port Moresby	Papua New Guinea	POM	AYPY	-9.443380	147.220001	146	10	U	Pacific/Port_Moresby	airport	OurAirports

# Create Point geometries from airport longitude and latitude
# The zip() function in Python combines multiple iterables (like lists or columns) into tuples, pairing elements by their positions.
airport_geometry = [Point(xy) for xy in zip(airports['long'],
                                            airports['lat'])]

# Create a GeoDataFrame from the airports data
airport_geodata = gpd.GeoDataFrame(airports,
                                   crs="EPSG:4326",
                                   geometry=airport_geometry)

# Create a quick plot of our data
fig, ax = plt.subplots(facecolor='black',
                       subplot_kw={'projection': ccrs.Robinson()},
                       figsize=(20,20))
ax.patch.set_facecolor('black')

airport_geodata.plot(ax=ax, transform=ccrs.PlateCarree(),
                     markersize=4, alpha=1, color='crimson',
                     edgecolors='none')
ax.set_ylim(-7000000, 9000000)
plt.show()

🌍 GeoPandas vs. Cartopy: When to Use Each for Map Visualization

GeoPandas and Cartopy are both powerful geospatial tools in Python, but they serve different purposes:

✅ GeoPandas + Matplotlib

• Best for local or regional maps

• Uses flat projections (e.g., UTM, Web Mercator)

• Quickly plots shapefiles and vector data

• Manual CRS transformation via .to_crs()

import geopandas as gpd

gdf = gdf.to_crs(epsg=3857)
gdf.plot()

🧭 Cartopy + Matplotlib

• Designed for global maps with curved projections (e.g., Robinson, Mollweide)

• Supports map features like coastlines, borders, land, and ocean

• Handles CRS transformations automatically

• Perfect for visualizing data on world-scale maps

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()

Use Case	GeoPandas	Cartopy
Flat, local maps	✅	✅
Global, curved projections	❌	✅
CRS transformations	✅ (manual)	✅ (automatic)
Coastlines & basemaps	❌	✅

🧭 Why Use projection and transform in Cartopy?

• ``projection=`` sets how the map is displayed, e.g., ccrs.Robinson().

• ``transform=`` tells Cartopy the CRS of your data, usually ccrs.PlateCarree() for lat/lon.

💡 Use both to ensure your geographic data appears correctly on the projected map.

Now let’s load the other data file containing flight information (flight routes).

routes = pd.read_csv(DATA_DIRECTORY / "routes.dat",
                     delimiter=',',
                     names=['airline', 'id', 'source_airport',
                            'source_airport_id', 'destination_airport',
                            'destination_airport_id', 'codeshare',
                            'stops', 'equitment'])
routes.head()

	airline	id	source_airport	source_airport_id	destination_airport	destination_airport_id	codeshare	stops	equitment
0	2B	410	AER	2965	KZN	2990	NaN	0	CR2
1	2B	410	ASF	2966	KZN	2990	NaN	0	CR2
2	2B	410	ASF	2966	MRV	2962	NaN	0	CR2
3	2B	410	CEK	2968	KZN	2990	NaN	0	CR2
4	2B	410	CEK	2968	OVB	4078	NaN	0	CR2

# Select relevant columns for source airports and create a new DataFrame
source_airports = airports[['name', 'iata', 'icao', 'lat', 'long']]

# Copy the source airports DataFrame to create a separate DataFrame for destination airports
destination_airports = source_airports.copy()

# Rename columns of the source airports DataFrame to indicate they are source attributes
source_airports.columns = [str(col) + '_source' for col in source_airports.columns]

# Rename columns of the destination airports DataFrame to indicate they are destination attributes
destination_airports.columns = [str(col) + '_destination' for col in destination_airports.columns]

# Filter the routes DataFrame to include only source and destination airport identifiers
routes = routes[['source_airport', 'destination_airport']]

# Merge the filtered routes DataFrame with source_airports on their IATA codes to add source airport details
routes = pd.merge(routes,
                  source_airports,
                  left_on='source_airport',  # Column in routes DataFrame to match
                  right_on='iata_source')    # Corresponding column in source_airports DataFrame

# Merge the updated routes DataFrame with destination_airports on their IATA codes to add destination airport details
routes = pd.merge(routes,
                  destination_airports,
                  left_on='destination_airport',  # Column in routes DataFrame to match
                  right_on='iata_destination')    # Corresponding column in destination_airports DataFrame

# Check the columns in our joined data
routes.columns

Index(['source_airport', 'destination_airport', 'name_source', 'iata_source',
       'icao_source', 'lat_source', 'long_source', 'name_destination',
       'iata_destination', 'icao_destination', 'lat_destination',
       'long_destination'],
      dtype='str')

Creating line objects to represent flight routes

Now that we did some data preparation, we use Shapely library to visualize flight routes between airports on a map. You can refresh your memory of Shapely library from AutoGIS course page Here’s a step-by-step breakdown of what each part of the code is doing:

1. Import LineString from Shapely: We import it to represent the flight routes geometrically.

2. Create LineString Objects for Each Route:

   • We iterate over each row in the routes DataFrame, which contains information about each flight route, including the source and destination airports’ coordinates (long_source, lat_source, long_destination, lat_destination).

   • For every route, we create a LineString object that connects the source airport to the destination airport using their longitude and latitude coordinates. This draws a line between the two points on the map.

3. Create a GeoDataFrame for Routes:

   • With the list of LineString objects (each representing a flight route), we create a new GeoDataFrame called routes_geodata. This GeoDataFrame contains all the original route information, plus a new geometry column where each entry is a LineString object corresponding to that route.

   • We also specify the coordinate reference system (CRS), which is ‘EPSG:4326’.

By transforming our route data into a GeoDataFrame with line geometries, we can easily visualize these flight paths on a map now.

from shapely.geometry import LineString

routes_geometry = [LineString([[routes.iloc[i]['long_source'],
                                routes.iloc[i]['lat_source']],
                               [routes.iloc[i]['long_destination'],
                                routes.iloc[i]['lat_destination']]])
                   for i in range(routes.shape[0])]

routes_geodata = gpd.GeoDataFrame(routes,
                                  geometry=routes_geometry,
                                  crs='EPSG:4326')

Now let’s make a quick plot to see what we have:

fig, ax = plt.subplots(figsize=(20,20))
ax.patch.set_facecolor('black')

routes_geodata.plot(ax=ax, color='white', linewidth=0.1)

plt.show()

Now here is where Cartopy gets to shine; using a geodetic coordinate system instead of the simple flat projection as above, can make our map nicer.

What is ccrs.Geodetic() in Cartopy?

ccrs.Geodetic() tells Cartopy that your data is in geodetic coordinates — standard latitude and longitude.

When used with transform=ccrs.Geodetic(), it ensures that lines (like flight routes) follow great-circle paths, the shortest path on a curved Earth.

Tip: Use it when plotting lat/lon lines on a projected map (e.g., Robinson) to get realistic curved routes instead of straight lines.

# Create a figure and axis with a Robinson projection
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.Robinson()},
                       figsize=(20,20))  # Set the size of the map

ax.patch.set_facecolor('black')  # Set the background color of the map to black

# Plot the routes on the map
routes_geodata.plot(ax=ax,
                    transform=ccrs.Geodetic(),  # Ensure the data uses a geodetic (lat/lon) projection
                    color='white',  # Set the color of the routes to white
                    linewidth=0.1,  # Set the width of the lines representing the routes
                    alpha=0.1)  # Set the transparency of the routes

# Set the y-axis limits to adjust the visible area of the map projection
ax.set_ylim(-7000000, 8800000)

# Display the plot
plt.show()

Let’s now add the aiport locations to the map. At the same time, we would like to show the number of flights per each airport. Why not do both at the same time using a proportional symbol map? For doing so, we need to do a bit more data preparation.

# Count the number of occurrences of each airport as a source airport
airport_source_count = routes.source_airport.value_counts()

# Count the number of occurrences of each airport as a destination airport
airport_destination_count = routes.destination_airport.value_counts()

# Display the counts of source airports
airport_source_count.head()

source_airport
ATL    915
ORD    558
LHR    527
PEK    525
CDG    524
Name: count, dtype: int64

# Convert the series of source airport counts into a DataFrame with two columns: 'airport' and 'source_airport_count'
airport_source_count = pd.DataFrame({'airport': airport_source_count.index,
                                     'source_airport_count': airport_source_count.values})

# Convert the series of destination airport counts into a DataFrame with two columns: 'airport' and 'destination_airport_count'
airport_destination_count = pd.DataFrame({'airport': airport_destination_count.index,
                                          'destination_airport_count': airport_destination_count.values})

# Merge the two DataFrames on the 'airport' column to combine source and destination counts for each airport
airport_counts = pd.merge(airport_source_count,
                          airport_destination_count,
                          on="airport")
# Calculate the total count (sum of source and destination counts) for each airport
airport_counts['count'] = airport_counts['source_airport_count'] + airport_counts['destination_airport_count']

# Merge the combined counts DataFrame with the airports DataFrame to add additional airport information
airport_counts = pd.merge(airport_counts,
                          airports,
                          left_on="airport",  # Key in airport_counts DataFrame
                          right_on="iata")    # Corresponding key in airports DataFrame

# Create Point geometries from airport longitude and latitude for geospatial visualization
geometry = [Point(xy) for xy in zip(airport_counts.long, airport_counts.lat)]

# Convert the airport_counts DataFrame into a GeoDataFrame
airport_counts = gpd.GeoDataFrame(airport_counts,
                                  geometry=geometry,
                                  crs="EPSG:4326")

# Calculate marker size based on the airport 'count' to visualize the airport traffic volume
airport_counts['markersize'] = airport_counts['count'] / 10  # Scale down the count for a suitable marker size

# Create a figure and axis with a Robinson projection, setting the size of the map
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.Robinson()},
                       figsize=(20,20))

# Set the background color of the map to black
ax.patch.set_facecolor('black')

# Plot the routes on the map with specific visual attributes
routes_geodata.plot(ax=ax,
                    transform=ccrs.Geodetic(),  # Ensure the data uses the Geodetic (latitude/longitude) projection
                    color='white',  # Set the routes color to white for contrast
                    linewidth=0.1,  # Specify the line width for the routes
                    alpha=0.1)  # Set the transparency level of the routes

# Plot the airports on the map, using the 'markersize' to represent traffic volume
airport_counts.plot(ax=ax,
                    transform=ccrs.PlateCarree(),  # Use PlateCarree projection for airport points
                    markersize=airport_counts['markersize'],  # Set marker size based on airport traffic volume
                    alpha=0.8,  # Set the transparency level of the airport markers
                    column=airport_counts['long'],  # Use the longitude for coloring (not typically meaningful, could be adjusted)
                    cmap='jet',  # Color map for the airports
                    edgecolors='none')  # No edge colors for the markers

# Adjust the y-axis limits to better fit the projection's display area
ax.set_ylim(-7000000, 8800000)

# Display the plot
plt.show()

Why Do transform=PlateCarree() and transform=Geodetic() Look Different?

Both refer to lat/lon (EPSG:4326), but:

• PlateCarree() assumes a flat grid, drawing straight lines between points.

• Geodetic() interprets points on a spherical Earth, drawing curved great-circle paths.

💡 transform= doesn’t transform your data — it tells Cartopy what CRS the data is in, so it can project it correctly onto the map.

Let’s still improve the map:

• We will add a legend to the map, enhancing interpretability by attempting to distinguish between routes and airports visually.

• Add a title to show to the viewer what the map is about.

• We are currently applying the colormap to the airport markers based on their longitude, adding a visual dimension to the data. (this could have also been used to still add one more variable to the map). Since this is not a crucial information in the map, let’s reduce its visual hierarchy by applying a more neutral color range coolwarm so it doesn’t dominate the visualization.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.lines import Line2D
from matplotlib.patches import Patch

# Initialize a figure with a Robinson projection and set the figure size and background color
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.Robinson()}, figsize=(20, 20))
ax.patch.set_facecolor('black')

# Plot the flight routes with specified visual attributes for a subtle background effect
routes_geodata.plot(ax=ax, transform=ccrs.Geodetic(), color='white', linewidth=0.1, alpha=0.1)

# Plot the airport counts with adjusted marker size and colormap
scatter = airport_counts.plot(ax=ax, transform=ccrs.PlateCarree(), markersize=airport_counts['markersize'],
                              alpha=0.8, column=airport_counts['long'], cmap='coolwarm', edgecolors='none')



# Adjust the vertical (y-axis) limits of the map to focus on the main areas of interest
ax.set_ylim(-7000000, 8800000)

# Manually create legend entries
legend_elements = [
    Line2D([0], [0], color='white', lw=2, alpha=0.6, label='Routes'),
    Line2D([0], [0], marker='o', color='w', markersize=10, markerfacecolor='grey', alpha=0.8, label='Airports', linestyle='None')
]

# Create the legend with the custom entries, specifying location and style without using 'title_color'
legend = ax.legend(handles=legend_elements, loc='upper left', facecolor='black', edgecolor='white', fontsize='large', title='Legend', title_fontsize='large')

# Set the legend title color after the legend has been created
legend.get_title().set_color("white")

# Set individual text colors for the legend items
legend.get_texts()[0].set_color('grey')  # Set color for the "Routes" legend item
legend.get_texts()[1].set_color('grey')  # Set color for the "Airports" legend item

# Add a title to the map
ax.set_title("Global Flight Routes and Airports", fontsize='xx-large', color='#db8a8a', loc='center', pad=10, fontweight='bold', family='sans-serif')

# save the map if you want
plt.savefig('Global_Flight_Routes.png', bbox_inches='tight', facecolor='white')

# Proceed to display the plot as before
plt.show()