Geospatial Visualization in Python

27. Geospatial Visualization in Python#

In this hands-on lecture you will put the principles from Communicating with Maps and Visualizations into code. You will:

See why rainbow colormaps mislead, and measure it.
Build choropleth maps that normalize, classify and project the data on purpose.
Make interactive maps with GeoPandas .explore() and Lonboard, and see why data size matters.
Look inside map tiles and PMTiles, read one over HTTP range requests, and write your own.
Fix a bad map.

All datasets are small files hosted on GitHub (Natural Earth via this repository, the Xarray tutorial data, and the PMTiles test files), so the notebook runs on a laptop.

Environment

Create the environment from this repository: run pixi install, then pixi run lab. It includes geopandas, matplotlib, xarray, mapclassify, cmcrameri, colorspacious, folium, lonboard, pmtiles, mapbox-vector-tile, mercantile, obstore and async-pmtiles. The Lonboard maps need a running Jupyter session; they do not show up on the course website.

You can access this notebook on this GitHub repo.

27.1. Setup#

import warnings
from pathlib import Path

warnings.filterwarnings("ignore", category=SyntaxWarning)  # an old docstring in colorspacious
warnings.filterwarnings("ignore", message="Numba not installed")  # mapclassify works without numba, just slower

import cmcrameri.cm as cmc
import geopandas as gpd
import mapclassify
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
import pooch
import shapely
import xarray as xr
from colorspacious import cspace_convert
from matplotlib import cbook

pooch.get_logger().setLevel("WARNING")
DATA_DIR = Path("data")
NE = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/"


def fetch(url, fname):
    """Download a file once and return its local path."""
    return pooch.retrieve(url, known_hash=None, fname=fname, path=DATA_DIR, progressbar=False)


countries_path = fetch(NE + "ne_110m_admin_0_countries.geojson", "ne_110m_admin_0_countries.geojson")
places_path = fetch(NE + "ne_10m_populated_places_simple.geojson", "ne_10m_populated_places_simple.geojson")
roads_path = fetch(NE + "ne_10m_roads.geojson", "ne_10m_roads.geojson")
air_path = fetch("https://raw.githubusercontent.com/pydata/xarray-data/master/air_temperature.nc", "air_temperature.nc")

27.3. Choropleth Maps#

27.3.1. Counts or rates?#

We use Natural Earth country boundaries with their population estimates. Antarctica is dropped. To compute areas and densities we first project to Equal Earth (EPSG:8857), an equal-area projection.

world = gpd.read_file(countries_path)
world = world.loc[world["NAME"] != "Antarctica", ["NAME", "ISO_A3", "CONTINENT", "POP_EST", "GDP_MD", "geometry"]]
world_ea = world.to_crs("EPSG:8857")
world_ea["area_km2"] = world_ea.area / 1e6
world_ea["density"] = world_ea["POP_EST"] / world_ea["area_km2"]
world_ea.sort_values("density", ascending=False).head()

	NAME	ISO_A3	CONTINENT	POP_EST	GDP_MD	geometry	area_km2	density
99	Bangladesh	BGD	Asia	163046161.0	302571	POLYGON ((8568892.999 2796633.753, 8586671.447...	133712.482996	1219.378755
79	Palestine	PSE	Asia	4685306.0	16276	POLYGON ((3151667.493 3941252.352, 3111811.269...	5036.412858	930.286323
140	Taiwan	TWN	Asia	23568378.0	1127000	POLYGON ((11169405.174 3086232.747, 11176611.6...	34311.270861	686.899010
77	Lebanon	LBN	Asia	6855713.0	51991	POLYGON ((3161493.926 4151743.847, 3138022.753...	10096.043450	679.049475
169	Rwanda	RWA	Africa	12626950.0	10354	POLYGON ((2913872.815 -145784.721, 2951564.335...	23364.503066	540.433065

fig, axes = plt.subplots(1, 2, figsize=(14, 4), layout="constrained")
world_ea.plot(ax=axes[0], column="POP_EST", cmap="YlGnBu", legend=True, legend_kwds={"label": "People", "shrink": 0.7})
axes[0].set_title("Total population")
world_ea.plot(ax=axes[1], column="density", cmap="YlGnBu", legend=True, legend_kwds={"label": "People per km²", "shrink": 0.7})
axes[1].set_title("Population density, unclassed")
for ax in axes:
    ax.set_axis_off()

../_images/2c149e5301b256fbf5aa368f1c183cb3a327ae271c75bd1ab5c652f012dc36fb.png

The left map says “China and India”, and almost nothing else. The density map is better, but it has a different problem: a handful of very dense countries (Bangladesh, small island and city states) stretch the color ramp, so almost every other country looks the same pale color. Look at the distribution:

ax = world_ea["density"].plot.hist(bins=60, color="0.6", figsize=(6, 3))
ax.set_xlabel("People per km²")
ax.set_ylabel("Number of countries");

../_images/6fd38c7f48bd93f47c8f6c6c4a0ea10e1f6d5dfbde016ae89a5ac59248256217.png

27.3.2. Classification#

With skewed data like this, a classed map usually works better. GeoPandas passes the scheme argument to mapclassify. Let’s compare three schemes with five classes each.

schemes = ["EqualInterval", "Quantiles", "FisherJenks"]
fig, axes = plt.subplots(1, 3, figsize=(16, 4.5), layout="constrained")
for ax, scheme in zip(axes, schemes):
    world_ea.plot(ax=ax, column="density", scheme=scheme, k=5, cmap="YlGnBu", edgecolor="white", linewidth=0.2,
                  legend=True, legend_kwds={"title": "People per km²", "fmt": "{:.0f}", "loc": "lower left", "fontsize": 7})
    ax.set_title(scheme)
    ax.set_axis_off()

../_images/5f6cc3272bb34780d2f5e1a9867be08d278d66b371a308d5071525ec781c0c95.png

for scheme in schemes:
    c = getattr(mapclassify, scheme)(world_ea["density"], k=5)
    print(f"{scheme:>14}: upper bounds {np.round(c.bins).astype(int)}, countries per class {c.counts}")

 EqualInterval: upper bounds [ 244  488  732  976 1219], countries per class [153  17   4   1   1]
     Quantiles: upper bounds [  19   51   89  167 1219], countries per class [36 35 35 35 35]
   FisherJenks: upper bounds [  71  193  395  687 1219], countries per class [89 59 19  7  2]

Question

Which scheme puts 153 of 176 countries in the lowest class? Which one would you use to show “where people are densely concentrated”, and which to rank countries into fifths? There is no single right answer, but you should be able to defend your choice and name it in the legend.

27.3.3. Projection#

Web Mercator (EPSG:3857) is the projection of almost every web map. Let’s compare areas measured in Web Mercator with true areas from Equal Earth.

world_wm = world.to_crs("EPSG:3857")
cmp = world_ea.set_index("NAME")[["area_km2"]].assign(area_webmercator_km2=world_wm.set_index("NAME").area / 1e6)
cmp["inflation"] = cmp["area_webmercator_km2"] / cmp["area_km2"]
cmp.loc[["Greenland", "Canada", "United States of America", "Brazil", "Dem. Rep. Congo", "Indonesia"]].round(1)

	area_km2	area_webmercator_km2	inflation
NAME
Greenland	2206878.9	36285503.0	16.4
Canada	10033994.2	52166480.7	5.2
United States of America	9513248.9	21862279.6	2.3
Brazil	8508179.1	9059541.4	1.1
Dem. Rep. Congo	2323379.2	2358239.7	1.0
Indonesia	1817750.9	1840122.0	1.0

fig, axes = plt.subplots(1, 2, figsize=(14, 4.5), layout="constrained")
for ax, gdf, title in [(axes[0], world_wm, "Web Mercator (EPSG:3857)"), (axes[1], world_ea, "Equal Earth (EPSG:8857)")]:
    gdf.plot(ax=ax, color="0.85", edgecolor="white", linewidth=0.3)
    gdf[gdf["NAME"].isin(["Greenland", "Dem. Rep. Congo"])].plot(ax=ax, color=cmc.batlow(0.3))
    ax.set_title(title)
    ax.set_axis_off()

../_images/111d3fa7299e7d0adee04ee481db605f408bc591dcc1beb7b0e35520028046af.png

In Web Mercator, Greenland looks much larger than the Democratic Republic of the Congo. In reality the Congo is slightly larger. For choropleths, density maps, and any area calculation, use an equal-area CRS.

27.3.4. A finished map#

Putting it together: a rate, a stated classification, an equal-area projection, a title with the message, a legend with units, and the source.

fig, ax = plt.subplots(figsize=(11, 5.5))
world_ea.plot(ax=ax, column="density", scheme="FisherJenks", k=5, cmap="YlGnBu",
              edgecolor="white", linewidth=0.3, legend=True,
              legend_kwds={"title": "People per km²\n(natural breaks)", "fmt": "{:.0f}", "loc": "lower left", "frameon": False})
ax.set_title("Population is most concentrated in South Asia, East Asia and Western Europe", loc="left", fontsize=13)
ax.set_axis_off()
fig.text(0.13, 0.08, "Data: Natural Earth (1:110m), population estimates. Projection: Equal Earth (EPSG:8857).", fontsize=8, color="0.4");

../_images/44ac9bb00e36e9ddee2dd9a21b0b08d851ec51cf75318cea016a1b0dd2861658.png

27.4. Interactive Maps#

27.4.1. Small data: `.explore()`#

For small datasets, GeoDataFrame.explore() creates an interactive Leaflet map with folium. It accepts the same column, scheme and cmap arguments as .plot(). Hover over a country to see its attributes. Note that .explore() reprojects to Web Mercator for display, because that is what the basemap tiles use.

world_ea.explore(column="density", scheme="FisherJenks", k=5, cmap="YlGnBu",
                 tooltip=["NAME", "POP_EST", "density"], legend_kwds={"caption": "People per km²"})

Make this Notebook Trusted to load map: File -> Trust Notebook

27.4.2. When data gets bigger#

Folium converts the whole GeoDataFrame to GeoJSON text and sends it to the browser. That is fine for 176 countries. Let’s look at a larger dataset: the Natural Earth roads (1:10m).

# read only the columns we need
roads = gpd.read_file(roads_path, columns=["type", "name"])
places = gpd.read_file(places_path, columns=["name", "adm0name", "pop_max"])

n_coords = shapely.get_num_coordinates(roads.geometry).sum()
print(f"{len(roads):,} road features with {n_coords:,} vertices")

roads.to_parquet(DATA_DIR / "roads.parquet")
print(f"GeoJSON:    {Path(roads_path).stat().st_size / 1e6:6.1f} MB")
print(f"GeoParquet: {(DATA_DIR / 'roads.parquet').stat().st_size / 1e6:6.1f} MB")

56,600 road features with 709,166 vertices
GeoJSON:      50.5 MB
GeoParquet:    8.4 MB

About 57,000 lines and 700,000 vertices: 50 MB as GeoJSON text, but a fraction of that as GeoParquet. Calling roads.explore() would push all 50 MB of text into your browser and can freeze it (try it at your own risk). This is where Lonboard helps.

27.4.3. Lonboard#

The quickest way in is viz(), which chooses layer types for you. Let’s try it on the 7,300 populated places. Every Lonboard map keeps its own copy of its data, in Python and in your browser, so we will draw the roads only once, in the styled map below.

from lonboard import Map, PathLayer, ScatterplotLayer, viz
from lonboard.basemap import CartoStyle, MaplibreBasemap
from lonboard.colormap import apply_categorical_cmap, apply_continuous_cmap

viz(places)

For more control, build the layers yourself. Each styling property can be a single value or one value per feature, given as a NumPy array. Here we color roads by type with a qualitative (Okabe-Ito) palette, and draw cities as circles whose area is proportional to population (proportional symbols, since population is a count) and whose color shows the same value on a log scale with batlow.

road_colors = {
    "Major Highway": [213, 94, 0],
    "Secondary Highway": [0, 114, 178],
    "Road": [150, 150, 150],
}
road_type = roads["type"].where(roads["type"].isin(road_colors.keys()), "Road")

roads_layer = PathLayer.from_geopandas(
    roads,
    get_color=apply_categorical_cmap(road_type, road_colors),
    get_width=np.where(road_type == "Major Highway", 1.5, 0.7),
    width_units="pixels",
)

pop = places["pop_max"].clip(lower=1_000).to_numpy()
norm = mcolors.LogNorm(vmin=1e4, vmax=pop.max(), clip=True)
cities_layer = ScatterplotLayer.from_geopandas(
    places,
    get_fill_color=apply_continuous_cmap(norm(pop), cmc.batlow, alpha=0.85),
    get_radius=np.sqrt(pop) * 15,  # area proportional to population
    radius_units="meters",
    radius_min_pixels=1.5,
    stroked=False,
)

m = Map([roads_layer, cities_layer], basemap=MaplibreBasemap(style=CartoStyle.Positron), height=550)
m

Click a road or a city to see its attributes. Zoom in on Worcester and notice how smoothly it pans with 57,000 lines on screen.

Performance tip

Changing a property of an existing layer only sends that property to the browser. Re-creating the map sends all the data again. Run the next cell and watch the map above change.

cities_layer.get_fill_color = apply_continuous_cmap(norm(pop), cmc.lajolla, alpha=0.85)
roads_layer.opacity = 0.5

Finally, save the map as a standalone HTML file you can share or embed in a website:

m.to_html("roads_and_cities.html", title="Roads and populated places (Natural Earth)")

Memory

Each Lonboard map keeps a copy of its data in the Python kernel and in your browser tab, and it stays there as long as the map’s output is on the page. With large data, make one map and change its layers instead of making many maps, and close maps you no longer need with m.close() (or clear the cell output). Also keep JupyterLab’s Settings > Save Widget State Automatically turned off: otherwise all of that data is written into the .ipynb file when you save.

27.5. Tiles and PMTiles#

27.5.1. The tile pyramid#

The mercantile package does the tile arithmetic for the Web Mercator grid. Which tile contains Clark University at each zoom level, and how many tiles are there in the world?

import mercantile

clark = (-71.8227, 42.2510)  # lon, lat
for z in [0, 4, 8, 12, 16, 18]:
    t = mercantile.tile(*clark, z)
    print(f"z={z:>2}  tile {t.z}/{t.x}/{t.y}   tiles in the world: {4**z:>15,}")

z= 0  tile 0/0/0   tiles in the world:               1
z= 4  tile 4/4/5   tiles in the world:             256
z= 8  tile 8/76/94   tiles in the world:          65,536
z=12  tile 12/1230/1516   tiles in the world:      16,777,216
z=16  tile 16/19693/24266   tiles in the world:   4,294,967,296
z=18  tile 18/78772/97065   tiles in the world:  68,719,476,736

How many tiles would you need to cover Massachusetts at different zoom levels? This is roughly what a tile generator has to create.

ma_bbox = (-73.51, 41.24, -69.93, 42.89)
for z in [6, 10, 14, 16]:
    n = sum(1 for _ in mercantile.tiles(*ma_bbox, zooms=z))
    print(f"z={z:>2}: {n:>9,} tiles")

z= 6:         2 tiles
z=10:        96 tiles
z=14:    16,728 tiles
z=16:   265,118 tiles

27.5.2. Reading a PMTiles file with HTTP range requests#

The PMTiles project keeps a small vector tile archive of Florence, Italy, built from OpenStreetMap, in its GitHub repository. GitHub serves files with support for HTTP range requests, so we can treat it like a cloud bucket.

We give the pmtiles reader a function that fetches a byte range, and we count every byte we download.

import requests
from pmtiles.reader import Reader

URL = "https://raw.githubusercontent.com/protomaps/PMTiles/main/spec/v3/protomaps(vector)ODbL_firenze.pmtiles"
file_size = int(requests.head(URL, allow_redirects=True).headers["Content-Length"])

downloaded = []


def get_bytes(offset, length):
    """Fetch `length` bytes starting at `offset` with an HTTP range request."""
    r = requests.get(URL, headers={"Range": f"bytes={offset}-{offset + length - 1}"})
    r.raise_for_status()
    downloaded.append(length)
    return r.content


reader = Reader(get_bytes)
header = reader.header()
print(f"File size: {file_size / 1e6:.1f} MB")
print(f"Tile type: {header['tile_type'].name}, zoom {header['min_zoom']} to {header['max_zoom']}, "
      f"{header['addressed_tiles_count']} tiles")
print(f"Root directory starts at byte {header['root_offset']} (the header is 127 bytes)")
print(f"Layers: {[layer['id'] for layer in reader.metadata()['vector_layers']]}")

File size: 6.6 MB
Tile type: MVT, zoom 0 to 15, 92 tiles
Root directory starts at byte 127 (the header is 127 bytes)
Layers: ['boundaries', 'buildings', 'earth', 'landcover', 'landuse', 'natural', 'physical_line', 'physical_point', 'places', 'pois', 'roads', 'transit', 'water']

Now fetch the single zoom-14 tile that contains the Duomo of Florence:

z = 14
duomo = mercantile.tile(11.2560, 43.7731, z)
tile_bytes = reader.get(duomo.z, duomo.x, duomo.y)

print(f"Tile {duomo.z}/{duomo.x}/{duomo.y}: {len(tile_bytes) / 1e3:.0f} KB (compressed)")
print(f"{len(downloaded)} range requests, {sum(downloaded) / 1e3:.0f} KB downloaded "
      f"= {100 * sum(downloaded) / file_size:.1f}% of the file")

Tile 14/8704/5971: 84 KB (compressed)
6 range requests, 107 KB downloaded = 1.6% of the file

The reader fetched the header, the directory and one tile, and nothing else. A map in your browser does exactly the same thing, for the 10 to 20 tiles that cover the screen.

27.5.3. What is inside a vector tile?#

A vector tile is a small binary (Protocol Buffers) file in the Mapbox Vector Tile format, here also gzip compressed. Let’s decode it.

import gzip

import mapbox_vector_tile

tile = mapbox_vector_tile.decode(gzip.decompress(tile_bytes))
for layer_name, layer in tile.items():
    print(f"{layer_name:>15}: {len(layer['features']):4d} features, extent {layer['extent']}")

      buildings:   16 features, extent 4096
          earth:    1 features, extent 4096
        landuse:  201 features, extent 4096
        natural:    1 features, extent 4096
  physical_line:    1 features, extent 4096
 physical_point:    2 features, extent 4096
         places:    4 features, extent 4096
           pois:  226 features, extent 4096
          roads:  590 features, extent 4096
        transit:   83 features, extent 4096
          water:    5 features, extent 4096

Each layer has an extent of 4096: geometry coordinates are integers on a 4096 by 4096 grid inside the tile, not longitudes and latitudes. Let’s draw the roads and buildings layers of this one tile in tile coordinates.

from shapely.geometry import shape

fig, ax = plt.subplots(figsize=(6, 6))
for layer_name, color, kw in [("buildings", cmc.batlow(0.75), {"edgecolor": "none"}), ("roads", "0.2", {"linewidth": 0.6})]:
    geoms = [shape(f["geometry"]) for f in tile[layer_name]["features"]]
    gpd.GeoSeries(geoms).plot(ax=ax, color=color, **kw)
ax.set_xlim(0, 4096)
ax.set_ylim(0, 4096)
ax.set_title(f"Tile {duomo.z}/{duomo.x}/{duomo.y} in tile coordinates");

../_images/b5cf178e32089624d77a50a588f10007bc3ab77f6c6f56f060be99368c93c6f1.png

Notice that roads stop at the tile edge: features are clipped to the tile (plus a small buffer). The same road continues in the neighboring tile as a separate piece. At lower zoom levels, geometries are also simplified and small features dropped. That is why tiles are a format for display and not for analysis: you cannot compute a road’s length or a building’s area reliably from tiles.

27.5.4. Make your own PMTiles#

GDAL (3.8 and later) can write PMTiles, so GeoPandas can too. We write the roads for zoom levels 0 to 7. This takes about half a minute.

pmtiles_path = DATA_DIR / "roads.pmtiles"
pmtiles_path.unlink(missing_ok=True)  # the PMTiles driver cannot overwrite an existing file
roads.to_file(pmtiles_path, driver="PMTiles", layer="roads", MINZOOM=0, MAXZOOM=7)
print(f"roads.pmtiles: {pmtiles_path.stat().st_size / 1e6:.1f} MB")

roads.pmtiles: 10.7 MB

Open the local file with the reader, this time from disk, and count the features in the tile containing Worcester at different zoom levels.

from pmtiles.reader import MmapSource

with open(pmtiles_path, "rb") as f:
    local = Reader(MmapSource(f))
    for z in [2, 4, 7]:
        t = mercantile.tile(-71.80, 42.26, z)
        data = local.get(t.z, t.x, t.y)
        layer = mapbox_vector_tile.decode(gzip.decompress(data))["roads"]
        print(f"z={z}: tile {t.z}/{t.x}/{t.y} is {len(data) / 1e3:6.0f} KB with {len(layer['features']):5,} features")

z=2: tile 2/1/1 is    150 KB with 8,281 features
z=4: tile 4/4/5 is     40 KB with 1,786 features
z=7: tile 7/38/47 is      5 KB with   219 features

The zoom 2 tile covers a quarter of the Western Hemisphere (eastern North America, the Caribbean and most of South America) and holds thousands of road features. GDAL clips and simplifies geometries but keeps every feature at every zoom, so low-zoom tiles get large. Tools built for this job, such as tippecanoe, drop and merge features at low zoom levels to keep tiles small (a common target is under 500 KB per tile). For large datasets, use tippecanoe.

Try it

Go to pmtiles.io and drag your data/roads.pmtiles file onto the page. The viewer reads the file in your browser (nothing is uploaded) and lets you inspect the tiles and the metadata, zoom level by zoom level.

27.5.5. Display PMTiles on a map#

Vector PMTiles are drawn by MapLibre GL JS, the open-source web map library, with the small pmtiles plugin that turns the file into range requests. A complete web map is just a short HTML page: load the two libraries, register the pmtiles:// protocol, and give MapLibre a style that says which layers of the file to draw and how. The browser reads the remote file directly, and the style (colors, widths) is applied in the browser, so you can change it without rebuilding the tiles.

We save the page as firenze_map.html (open it in your browser: a static page plus a static PMTiles file is a web map with no server) and also show it below.

Why an iframe?

The page is shown inside its own frame (an iframe), isolated from the rest of the notebook. That matters here: Lonboard and leafmap’s MapLibre widget each bundle a different version of deck.gl, and deck.gl refuses to load twice in one page. In a notebook without Lonboard maps, leafmap.maplibregl.Map().add_pmtiles(...) does the same job as this cell in one line.

import json
from html import escape

from IPython.display import HTML

style = {
    "version": 8,
    "sources": {"firenze": {"type": "vector", "url": "pmtiles://" + URL,
                            "attribution": "© OpenStreetMap contributors"}},
    "layers": [
        {"id": "background", "type": "background", "paint": {"background-color": "#f2f2f2"}},
        {"id": "water", "source": "firenze", "source-layer": "water", "type": "fill",
         "paint": {"fill-color": "#a6cee3"}},
        {"id": "buildings", "source": "firenze", "source-layer": "buildings", "type": "fill",
         "paint": {"fill-color": "#e69f00", "fill-opacity": 0.6}},
        {"id": "roads", "source": "firenze", "source-layer": "roads", "type": "line",
         "paint": {"line-color": "#333333", "line-width": 1}},
    ],
}

page = f"""<!DOCTYPE html>
<html><head><meta charset="utf-8">
<link rel="stylesheet" href="https://unpkg.com/maplibre-gl@5.24.0/dist/maplibre-gl.css">
<script src="https://unpkg.com/maplibre-gl@5.24.0/dist/maplibre-gl.js"></script>
<script src="https://unpkg.com/pmtiles@4.5.0/dist/pmtiles.js"></script>
<style>html, body, #map {{margin: 0; height: 100%;}}</style>
</head><body><div id="map"></div><script>
const protocol = new pmtiles.Protocol();
maplibregl.addProtocol("pmtiles", protocol.tile);
const map = new maplibregl.Map({{container: "map", style: {json.dumps(style)}, center: [11.256, 43.773], zoom: 14}});
map.addControl(new maplibregl.NavigationControl());
</script></body></html>"""

# A standalone web page: open it in any browser
Path("firenze_map.html").write_text(page)

# Show the same page here, inside its own frame
HTML(f'<div><iframe srcdoc="{escape(page)}" style="width: 100%; height: 450px; border: none;"></iframe></div>')

Raster PMTiles (images) can be displayed with Lonboard’s RasterLayer. The archive below holds USGS imagery of Mount Whitney, California, as WebP tiles for zoom levels 8 to 15. In this case Python reads the tiles with async-pmtiles and obstore and passes them to the map. Jupyter supports await at the top level of a cell.

from async_pmtiles import PMTilesReader
from lonboard import RasterLayer
from obstore.store import HTTPStore

store = HTTPStore("https://raw.githubusercontent.com/protomaps/PMTiles/main/spec/v3")
whitney = await PMTilesReader.open("usgs-mt-whitney-8-15-webp-512.pmtiles", store=store)
print(whitney.tile_type, whitney.minzoom, whitney.maxzoom)

Map(RasterLayer.from_pmtiles(whitney), height=450)

TileType.WEBP 8 15

27.6. Exercise: Fix This Map#

A colleague sends you the map below to show “how wealthy countries are”. List everything wrong with it, then make a better version in the empty cell. Think about: counts versus rates, the colormap, classification, projection, labels, and the source.

ax = world.to_crs("EPSG:3857").plot(column="GDP_MD", cmap="jet", legend=True, figsize=(10, 5))
ax.set_title("gdp_md");

../_images/b4dfef1683538fcfe951d2cf330361126fffe0665b778aae1e0d4e41ae23aa0f.png

# Your improved map here

Geospatial Visualization in Python

Contents

27. Geospatial Visualization in Python#

27.1. Setup#

27.2. Color#

27.2.1. The same data, three colormaps#

27.2.2. Measure it: perceived lightness#

27.2.3. Color vision deficiency#

27.2.4. Diverging data: center the colormap#

27.2.5. Comparing panels: use one scale#

27.2.6. The same rules apply to charts#

27.3. Choropleth Maps#

27.3.1. Counts or rates?#

27.3.2. Classification#

27.3.3. Projection#

27.3.4. A finished map#

27.4. Interactive Maps#

27.4.1. Small data: `.explore()`#

27.4.2. When data gets bigger#

27.4.3. Lonboard#

27.5. Tiles and PMTiles#

27.5.1. The tile pyramid#

27.5.2. Reading a PMTiles file with HTTP range requests#

27.5.3. What is inside a vector tile?#

27.5.4. Make your own PMTiles#

27.5.5. Display PMTiles on a map#

27.6. Exercise: Fix This Map#

Geospatial Visualization in Python

Contents

27. Geospatial Visualization in Python#

27.1. Setup#

27.2. Color#

27.2.1. The same data, three colormaps#

27.2.2. Measure it: perceived lightness#

27.2.3. Color vision deficiency#

27.2.4. Diverging data: center the colormap#

27.2.5. Comparing panels: use one scale#

27.2.6. The same rules apply to charts#

27.3. Choropleth Maps#

27.3.1. Counts or rates?#

27.3.2. Classification#

27.3.3. Projection#

27.3.4. A finished map#

27.4. Interactive Maps#

27.4.1. Small data: .explore()#

27.4.2. When data gets bigger#

27.4.3. Lonboard#

27.5. Tiles and PMTiles#

27.5.1. The tile pyramid#

27.5.2. Reading a PMTiles file with HTTP range requests#

27.5.3. What is inside a vector tile?#

27.5.4. Make your own PMTiles#

27.5.5. Display PMTiles on a map#

27.6. Exercise: Fix This Map#

27.4.1. Small data: `.explore()`#