Skip to main content
Geographical Plots

Geographical Plots #

Warning: This post hasn't been updated for over a year. The information may be out of date.

(Maybe) Useful resources #

Courses:

Things to check out later:
Things to not check out later:

Basic plotting flow #

  1. Have data (to plot) with latitude and longitude columns;
  2. Get map (vector map/shape file/polygons/tiles/webmap) of the desired range; (some packages, like Basemap, seem to come with global and US maps, but not of other countries)
  3. Import both datasets;
  4. Plot map data first, then layer the actual data on top of it.

Main usage of some packages #

I will be mainly using matplotlib and seaborn to plot with cartopy. If I got extra time, I will try out geoplot as well.

geopandas
Pandas DataFrame equivalence for data or map that are geographical polygons. Also plots the polygons and points.
Introduction to GeoPandas — GeoPandas 0+untagged.50.g5558c35.dirty documentation
contextily
Auto-add webmap to data plot as background image.
Introduction guide to contextily — contextily 1.1.0 documentation
xyzservices
API interfaces that contextily uses. Sometimes debugging contextily unavoidably leads to it, for example, looking for the now gone Stamen provider.
Gallery - xyzservices 2023.10.1
cartopy
Handle map projection of plots (Matplotlib axes) and data to be plotted on those axes.
Gallery — cartopy 0.22.0 documentation
matplotlib/seaborn
My go-to Python plotting packages. See Matplotlib/Seaborn.
geoplot
Claims to be a geospatial version of seaborn. May try later.
Gallery — geoplot 0.5.0 documentation
Others
I was told that folium, plotly, bokeh, geoviews and many other plotting packages could plot geospatial figures as well. But since I am already familiar with matplotlib and seaborn, I decided against trying other packages. I remember trying bokeh a few years ago, but did not stick with it for some reason.

How to get dataset with latitude/longitude #

Datasets with lat/lon:

Convert postcode to lat/long:

  • Worldwide by-country lookup TSVs: GeoNames (column names: ["country code", "postal code", "place name", "admin name1", "admin code1", "admin name2", "admin code2", "admin name3", "admin code3", "latitude", "longitude", "accuracy"], see readme.txt)
    • Python wrapper: symerio/pgeocode (does not support the “full” version for CA,NL and UK yet)
  • geopandas.tools.geocode: geopy/geopy wrapper. Most decent providers (like Google) require API keys.

How to get map #

Two methods, bring your own map (as polygons or image) and plot them manually, or use some packages to download and plot automatically.

I am a lazy person. So far, I’ve found two different automatic approach, using cartopy and contextily, respectively.

Adding map with cartopy #

It is not maps per se, but I will be abusing the “feature” function in cartopy as makeshift maps.

Refs:

Toy example with no projection:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import cartopy.feature as cfeature

# define extra features
states_provinces = cfeature.NaturalEarthFeature(
    category='cultural',
    name='admin_1_states_provinces_lines',
    scale='50m',
)
urban_areas = cfeature.NaturalEarthFeature(
    category='cultural',
    name='urban_areas',
    scale='50m',
)

# begin plotting
plt.clf()
fig,ax = plt.subplots()
sns.scatterplot(ax=ax, data=df, x="longitude", y="latitude", hue="Rating")

# start cartopy settings
ax.set_global()

ax.add_feature(cfeature.LAND, facecolor="silver")
# ax.add_feature(cfeature.OCEAN)
# ax.add_feature(cfeature.COASTLINE, edgecolor="silver")
ax.add_feature(cfeature.BORDERS, edgecolor="white")
# ax.add_feature(states_provinces, linestyle=':', edgecolor="white", facecolor="none")
# ax.add_feature(cfeature.LAKES)
# ax.add_feature(cfeature.RIVERS, edgecolor="white")
# ax.add_feature(urban_areas, facecolor="silver, alpha=0.5)

ax.gridlines(draw_labels=True, alpha=0.3)

# matplotlib settings
ax.legend(frameon=True)
plt.show()
ax.add_image

cartopy has its own ax.add_image and cartopy.io.img_tiles methods, but the documentation is very ambiguous and the method not easy to use.

Docs:

Toy example with no projection: (remember Stamen is no longer available)

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import cartopy.io.img_tiles as cimgt

# begin plotting
plt.clf()
fig,ax = plt.subplots()
sns.scatterplot(ax=ax, data=df, x="longitude", y="latitude", hue="Rating")

# cartopy settings
ax.set_global()
ax.add_image(cimgt.OSM(), 5) # the larger zoom level, the smaller the label fonts
ax.gridlines(draw_labels=True, alpha=0.3)

# matplotlib settings
ax.legend(frameon=True)
plt.show()

Adding map with contextily #

The maps are beautiful out of the box, but handling projection becomes a problem if you are not using geopandas.

Ref: Creating a Map with XYZ Tiles using Geopandas, Matplotlib, Contextily, and XYZServices | Tutorials/Post - Remote Sensing, GIS, Earth System, Geo-AI/ML

Toy example with cartopy and seaborn with no projection, which translates to equirectangular projection on plot axes, aka PlateCarree:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import pandas as pd
import contextily as ctx
import matplotlib.pyplot as plt
import seaborn as sns

# begin plotting
plt.clf()
fig,ax = plt.subplots()
sns.scatterplot(ax=ax, data=df, x="longitude", y="latitude", hue="Rating")

# add map to background
ctx.add_basemap(
    ax,
    source=ctx.providers.CartoDB.Positron,
    zoom="auto",
    crs="EPSG:4326",
)

# matplotlib settings
ax.legend()
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.show()

A more involved way of using the add_basemap function that I have yet to try: Adding a background map to plots — GeoPandas 0+untagged.50.g5558c35.dirty documentation

Arbitrary map projection with cartopy and seaborn #

Most examples online read as if only FacetGrid could use custom projections, which is not true. I think any plot could use any projection with a proper amount of matplotlib setup.

Refs:

projection vs transform #

Refs:

Example with cartopy features #

TBA

Calculations #

Distance #

Example from: python - Getting distance between two points based on latitude/longitude - Stack Overflow, my answer

import cartopy.geodesic as cgeo

coords_1 = (52.2296756, 21.0122287)
coords_2 = (52.406374, 16.9251681)
coords_3 = (52.406374, 10)

# Note: cartopy only accepts lon-lat pairs!!
globe = cgeo.Geodesic()

inv = globe.inverse(coords_1[::-1], coords_2[::-1])
print(inv.T[0]/1000)
# [279.3529016]

inv_2 = globe.inverse(coords_1[::-1], [coords_2[::-1], coords_3[::-1]])
print(inv_2.T[0]/1000)
# [279.3529016  750.45799898]

Wrapper function for easier usage with DataFrame:

import cartopy.geodesic as cgeo

def distance(longitudes, latitudes):
    assert(1 < len(longitudes))
    assert(1 < len(latitudes))
    assert(len(longitudes) == len(latitudes))

    places = list(zip(longitudes, latitudes))
    globe = cgeo.Geodesic()
    inv = globe.inverse(places[0], places[1:])
    return (inv.T[0]/1000)

Entertainment #

xkcd: Map Projections
977: Map Projections - explain xkcd
xkcd’s “Map Projections”, animated - YouTube gives an idea how much the lands are distorted when you spin the globe with each projection.

2025-04-06