Creating Maps

Data 304: Visualizing Data and Models

Airport data

R/altair
Python/Altair

altair::import_vega_data() can be use in R. Notice that URL can be recovered if you ever need to access the data that way.

library(altair)
library(vegabrite)

vega_data <-
  import_vega_data()

airports <- vega_data$airports()       # note the parens!
head(airports, 3) |> pander::pander()

Table continues below
iata	name	city	state	country	latitude
00M	Thigpen	Bay Springs	MS	USA	31.95
00R	Livingston Municipal	Livingston	TX	USA	30.69
00V	Meadow Lake	Colorado Springs	CO	USA	38.95

longitude
-89.23
-95.02
-104.6

vega_data$airports$description         # description of the data

[1] "This dataset lists US airports, including airport code, city, state, latitude, and longitude. This dataset is a subset of the data compiled and published at http://ourairports.com/data/, and is in the public domain."

vega_data$airports$url                 # URL if you want that

[1] "https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/airports.csv"

vega_data$airports$filepath            # Where the data live on your machine

[1] "/Users/rpruim/.virtualenvs/r-reticulate/lib/python3.12/site-packages/vega_datasets/_data/airports.csv"

vega_data$list_datasets()              # A list of all datasets. Use _ in place of - to access

 [1] "7zip"                            "airports"                       
 [3] "annual-precip"                   "anscombe"                       
 [5] "barley"                          "birdstrikes"                    
 [7] "budget"                          "budgets"                        
 [9] "burtin"                          "cars"                           
[11] "climate"                         "co2-concentration"              
[13] "countries"                       "crimea"                         
[15] "disasters"                       "driving"                        
[17] "earthquakes"                     "ffox"                           
[19] "flare"                           "flare-dependencies"             
[21] "flights-10k"                     "flights-200k"                   
[23] "flights-20k"                     "flights-2k"                     
[25] "flights-3m"                      "flights-5k"                     
[27] "flights-airport"                 "gapminder"                      
[29] "gapminder-health-income"         "gimp"                           
[31] "github"                          "graticule"                      
[33] "income"                          "iowa-electricity"               
[35] "iris"                            "jobs"                           
[37] "la-riots"                        "londonBoroughs"                 
[39] "londonCentroids"                 "londonTubeLines"                
[41] "lookup_groups"                   "lookup_people"                  
[43] "miserables"                      "monarchs"                       
[45] "movies"                          "normal-2d"                      
[47] "obesity"                         "ohlc"                           
[49] "points"                          "population"                     
[51] "population_engineers_hurricanes" "seattle-temps"                  
[53] "seattle-weather"                 "sf-temps"                       
[55] "sp500"                           "stocks"                         
[57] "udistrict"                       "unemployment"                   
[59] "unemployment-across-industries"  "uniform-2d"                     
[61] "us-10m"                          "us-employment"                  
[63] "us-state-capitals"               "volcano"                        
[65] "weather"                         "weball26"                       
[67] "wheat"                           "windvectors"                    
[69] "world-110m"                      "zipcodes"

The vega_datasets package in Python can simplify access to the vega datasets. Notice that URL can be recovered if you ever need to access the data that way.

import altair as alt
from vega_datasets import data
airports = data.airports()
print(data.airports.url)

https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/airports.csv

print(data.airports.description)

This dataset lists US airports, including airport code, city, state, latitude, and longitude. This dataset is a subset of the data compiled and published at http://ourairports.com/data/, and is in the public domain.

airports.head()

  iata                  name              city  ... country   latitude   longitude
0  00M               Thigpen       Bay Springs  ...     USA  31.953765  -89.234505
1  00R  Livingston Municipal        Livingston  ...     USA  30.685861  -95.017928
2  00V           Meadow Lake  Colorado Springs  ...     USA  38.945749 -104.569893
3  01G          Perry-Warsaw             Perry  ...     USA  42.741347  -78.052081
4  01J      Hilliard Airpark          Hilliard  ...     USA  30.688012  -81.905944

[5 rows x 7 columns]

Maps require projection

Need to “flatten the sphere”.

\[\mbox{projection} : \mbox{(lon, lat)} \to (x, y)\]

Many different projections

All projections distort something.
Some projections are designed to preserve a particular. feature (area, direction, distance, etc.).
Others make compromises between various trade-offs.
Some are designed for specific use cases.

Available projections

Vega-Lite uses projections from d3-geo library.

Some common projections

Type	Description
albers	The Albers’ equal-area conic projection. This is a U.S.-centric configuration of “conicEqualArea”.
albersUsa	A U.S.-centric composite with projections for the lower 48 states, Hawaii, and Alaska (scaled to 0.35 times the true relative area).
azimuthalEqualArea	The azimuthal equal-area projection.
azimuthalEquidistant	The azimuthal equidistant projection.
conicConformal	The conic conformal projection. The parallels default to [30°, 30°] resulting in flat top.
conicEqualArea	The Albers’ equal-area conic projection.
conicEquidistant	The conic equidistant projection.
equalEarth	The Equal Earth projection, by Bojan Šavrič et al., 2018.
equirectangular	The equirectangular (plate carrée) projection, akin to use longitude, latitude directly.
gnomonic	The gnomonic projection.
identity	The identity projection. Also supports additional boolean reflectX and reflectY parameters.
mercator	The spherical Mercator projection. Uses a default clipExtent such that the world is projected to a square, clipped to approximately ±85° latitude.
orthographic	The orthographic projection.
stereographic	The stereographic projection.
transverseMercator	The transverse spherical Mercator projection. Uses a default clipExtent such that the world is projected to a square, clipped to approximately ±85° latitude.

Code

world_map_url <- vega_data$world_110m$url

some_projections <- 
  c("albers",
    "albersUsa",
    "azimuthalEqualArea",
    "azimuthalEquidistant",
    "conicEqualArea",
    "conicEquidistant",
    "equalEarth",
    "equirectangular",
    "gnomonic",
    "mercator",
    "naturalEarth1",
    "orthographic",
    "stereographic",
    "transverseMercator")

vl_chart(width = 500, height = 300) |>
  vl_add_parameter("projection") |>
  vl_bind_select_input("projection", options = some_projections) |> 
  vl_add_data(
    url = world_map_url, 
    format = list(type = "topojson", feature = "countries")) |>
  vl_add_properties(projection = list(type = list(expr = "projection"))) |>
  vl_mark_geoshape(fill = "lightgray", stroke = "gray")

Code

import altair as alt
from vega_datasets import data

source = alt.topo_feature(data.world_110m.url, 'countries')

input_dropdown = alt.binding_select(options=[
    "albers",
    "albersUsa",
    "azimuthalEqualArea",
    "azimuthalEquidistant",
    "conicEqualArea",
    "conicEquidistant",
    "equalEarth",
    "equirectangular",
    "gnomonic",
    "mercator",
    "naturalEarth1",
    "orthographic",
    "stereographic",
    "transverseMercator"
], name='Projection ')
param_projection = alt.param(value="equalEarth", bind=input_dropdown)

(alt.Chart(source, width=500, height=300).mark_geoshape(
    fill='lightgray',
    stroke='gray'
).project(
    type=alt.expr(param_projection.name)
).add_params(param_projection)
)

Specifying a projection

vegabrite
Python/Altair
Vega-Lite

vl_chart() |>
  vl_add_properties(projection = list(type = "albers"))

Use .project() in altair.

alt.Chart().project(type = "albers")

The projection is specified in top level of JSON:

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "projection": {
    "type": "albers"
  }
}

Can also be configured in the “config” section of the specification.

// Top-level View Specification
{ ...,
  "config": {          // Configuration Object
    "projection": { "type": ..., ... },   // - Projection Configuration
    ...
  }
}

Map of airports

vegabrite
Pthon/Altair
Vega-Lite

vl_chart(width = 800, height = 500) |>
  vl_mark_circle(opacity = 0.5, size = 5) |>
  vl_encode(longitude = "longitude:Q", latitude = "latitude:Q") |>
  vl_add_properties(project = list(type = "mercator")) |>
  vl_add_data_url(vega_data$airports$url)

airports_chart = alt.Chart(airports).mark_circle(opacity = 0.5).encode(
    longitude='longitude:Q', latitude='latitude:Q',
    size=alt.value(10), tooltip='name'
).project("mercator")

airports_chart.properties( width=800, height=500 )

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "height": 500,
  "width": 800,
  "mark": {
    "opacity": 0.5,
    "size": 5,
    "type": "circle"
  },
  "encoding": {
    "longitude": {
      "field": "longitude",
      "type": "quantitative"
    },
    "latitude": {
      "field": "latitude",
      "type": "quantitative"
    }
  },
  "projection": {
    "type": "mercator"
  },
  "data": {
    "url": "https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/airports.csv"
  }
}

AlbersUsa projection

This unusual projection moves and rescales Alaska and Hawaii.

vegabrite
Pyhton/Altair
Vega-Lite

Code

airports_points <-
  vl_chart(width = 450, height = 300) |>
  vl_mark_circle(opacity = 0.5, size = 5) |>
  vl_encode(longitude = "longitude:Q", latitude = "latitude:Q") |>
  vl_add_properties(project = list(type = "albersUsa")) |>
  vl_add_data_url(vega_data$airports$url)
airports_points

Code

airports_chart = alt.Chart(airports).mark_circle(opacity = 0.5).encode(
    longitude='longitude:Q', latitude='latitude:Q',
    size=alt.value(10), tooltip='name'
).project("albersUsa")

airports_chart.properties( width=450, height=300 )

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "height": 300,
  "width": 450,
  "mark": {
    "opacity": 0.5,
    "size": 5,
    "type": "circle"
  },
  "encoding": {
    "longitude": {
      "field": "longitude",
      "type": "quantitative"
    },
    "latitude": {
      "field": "latitude",
      "type": "quantitative"
    }
  },
  "projection": {
    "type": "albersUsa"
  },
  "data": {
    "url": "https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/airports.csv"
  }
}

Geospatial data

High level description

collection of features
each feature is a point, line (e.g. Grand River), polygon (Kansas), or set of these (Michigan).
- list of coordinates given for each in lon/lat
features may have additional data properties as well.

Vega-Lite prefers GeoJSON: or TopoJSON: formats.
Other formats (like shapefiles) can be converted to one of these formats.

Map of US States

What is unusual (for Vega-Lite) about the code for this map?

vegabrite
Python/Altair
Vega-Lite

us_map_url <- vega_data$us_10m$url
state_map <-
  vl_chart() |>
  vl_add_data(
    url = us_map_url, 
    format = list(type = "topojson", feature = "states")) |>
  vl_add_properties(projection = list(type = "albersUsa")) |>
  vl_mark_geoshape(fill = "transparent", stroke = "steelblue") 
state_map |> vl_add_properties(width = 500, height = 300)

print(data.us_10m.url)

https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/us-10m.json

states = alt.topo_feature(data.us_10m.url, feature = 'states')

state_map = alt.Chart(states).mark_geoshape(
    fill = 'transparent',
    stroke = 'steelblue'
).project('albersUsa')

state_map.properties(width = 500, height = 300)

state_map |> format()

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "data": {
    "format": {
      "type": "topojson",
      "feature": "states"
    },
    "url": "https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/us-10m.json"
  },
  "projection": {
    "type": "albersUsa"
  },
  "mark": {
    "fill": "transparent",
    "stroke": "steelblue",
    "type": "geoshape"
  }
}

What is us-10m?

https://cdn.jsdelivr.net/npm/vega-datasets@2.11.0/data/us-10m.json

A (topoJSON) JSON object with geo data.

the states feature contains information describing the state boundaries

GeoJSON and TopoJSON

TopoJSON:

GeoJSON:

https://cdn.jsdelivr.net/npm/vega-datasets@1.29.0/data/earthquakes.json

Combining layers

vegabrite
Python/Altair
Vega-Lite

# Use + or vl_layer() to make a 2-layer plot
state_map + airports_points

(state_map + airports_chart).properties(width = 500, height = 400)

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "layer": [
    {
      "data": {
        "format": {
          "type": "topojson",
          "feature": "states"
        },
        "url": "https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/us-10m.json"
      },
      "projection": {
        "type": "albersUsa"
      },
      "mark": {
        "fill": "transparent",
        "stroke": "steelblue",
        "type": "geoshape"
      }
    },
    {
      "height": 300,
      "width": 450,
      "mark": {
        "opacity": 0.5,
        "size": 5,
        "type": "circle"
      },
      "encoding": {
        "longitude": {
          "field": "longitude",
          "type": "quantitative"
        },
        "latitude": {
          "field": "latitude",
          "type": "quantitative"
        }
      },
      "projection": {
        "type": "albersUsa"
      },
      "data": {
        "url": "https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/airports.csv"
      }
    }
  ]
}

County level data

TopoJSON and GeoJSON files can contain data at different levels (features). This time we will use the county borders rather than the state borders.

vegabrite
Python/Altair

us_map_url <- vega_data$us_10m$url
county_map <-
  vl_chart(width = 500, height = 300) |>
  vl_add_data(
    url = us_map_url, 
    format = list(type = "topojson", feature = "counties")) |>
  vl_add_properties(projection = list(type = "albersUsa")) |>
  vl_mark_geoshape(fill = "transparent", stroke = "lightgray") 

county_map

counties = alt.topo_feature(data.us_10m.url, feature='counties')
county_map = alt.Chart(counties).mark_geoshape(
    fill='#eeeeee',
    stroke='lightgray'
).project('albersUsa')

county_map.properties(width = 600, height = 400)

Layering maps (States and counties)

vegabrite
Python/altair

county_map + state_map

(county_map + state_map).properties(width = 600, height = 400)

How do we make a map like this?

Need data with unemployment figures by county
Need to combine this with county map data (lookup)
Need a way to identify the same county in the two data sets (data consistency)

Unemployment data

vegabrite
Python/Altair

vega_data$unemployment() |> head(3) |> pander::pander()

id	rate
1001	0.097
1003	0.091
1005	0.134

print(data.unemployment.url)

https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/unemployment.tsv

print(data.unemployment().head(3))

     id   rate
0  1001  0.097
1  1003  0.091
2  1005  0.134

Important: id is used in the same way as in the geo data

Choropleth maps with lookup

Remember: Lookup = left (outer) join.

Note: In this case ‘id’ is used in both the geospatial data and the auxiliary data, but in other cases the names might differ.

vegabrite
Python/Altair

county_map |>
  vl_lookup(
    lookup = "id",                  # name of key variable in "left" data
    from = 
      list(
        data = list(url = vega_data$unemployment$url), 
        key = "id",                 # name of key variable in "right" data
        fields = list("rate"))) |>  # variables to merge into left
  vl_encode_stroke(value = "transparent") |>
  vl_encode_fill("rate:Q") |>
  vl_legend_fill(title = "unemployment rate", format = ".0%", type = "gradient") |>
  vl_add_properties(width = 350, height = 200)

county_map.transform_lookup(
  lookup='id',
  from_=alt.LookupData(data.unemployment(), 'id', ['rate'])
  ).encode(
    fill = "rate:Q",
    stroke = "rate:Q"
    ).properties(width = 600, height = 400)

Inspecting TopoJSON/GeoJSON files

The world map data in the Vega-Lite data sets doesn’t include country names!

R
Python

world_geo <- vega_data$world_110m()
names(world_geo)

[1] "type"      "transform" "objects"   "arcs"

names(world_geo$objects)

[1] "land"      "countries"

names(world_geo$objects$countries)

[1] "type"       "geometries"

world_geo$objects$countries$geometries[[1]]

$type
[1] "Polygon"

$arcs
$arcs[[1]]
[1] 499 500 501 502 503 504


$id
[1] 4

world1 = data.world_110m()
print(type(world1))

<class 'dict'>

print(world1.keys())

dict_keys(['type', 'transform', 'objects', 'arcs'])

print(world1['objects'].keys())

dict_keys(['land', 'countries'])

print(world1['objects']['countries'].keys())

dict_keys(['type', 'geometries'])

print(world1['objects']['countries']['geometries'][0])

{'type': 'Polygon', 'arcs': [[499, 500, 501, 502, 503, 504]], 'id': 4}

Another source of map data:

github.com/topojson

R
Python

library(jsonlite)

world2_url <- 'https://cdn.jsdelivr.net/npm/world-atlas@2/countries-110m.json'
world2_geo <- fromJSON(world2_url)
names(world2_geo)

[1] "type"      "objects"   "arcs"      "bbox"      "transform"

names(world2_geo$objects$countries$geometries)

[1] "type"       "arcs"       "id"         "properties"

names(world2_geo$objects$countries$geometries$properties)

[1] "name"

world2_geo$objects$countries$geometries$properties$name[1:3]

[1] "Fiji"      "Tanzania"  "W. Sahara"

import json 
from urllib.request import urlopen

world2_url = 'https://cdn.jsdelivr.net/npm/world-atlas@2/countries-110m.json'

world2 = json.load(urlopen(world2_url))
print(type(world2))

<class 'dict'>

print(world2.keys())

dict_keys(['type', 'objects', 'arcs', 'bbox', 'transform'])

print(world2['objects']['countries'].keys())

dict_keys(['type', 'geometries'])

print(world2['objects']['countries']['geometries'][0])

{'type': 'MultiPolygon', 'arcs': [[[0]], [[1]]], 'id': '242', 'properties': {'name': 'Fiji'}}

Ah, there is the country name – in objects > countries > properties > names. We can access this as the field properties.name in Vega-Lite.

Comparing country names

Want to see all the country names? Here is a way:

R

country_names <- world2_geo$objects$countries$geometries$properties$name
head(country_names)

[1] "Fiji"                     "Tanzania"                
[3] "W. Sahara"                "Canada"                  
[5] "United States of America" "Kazakhstan"

Now let’s check if the names match what is in the Gapminder data set.

gapminder <- vega_data$gapminder()
common_names <- intersect(country_names, gapminder$country)
world2_only_names <- setdiff(country_names, gapminder$country)
gapminder_only_names <- setdiff(gapminder$country, country_names)

tibble::tibble(common = length(common_names), 
  world2_only = length(world2_only_names), 
  gapminder_only = length(gapminder_only_names)) |> pander::pander()

common	world2_only	gapminder_only
57	120	6

Python

country_names = [ p['properties']['name'] 
  for p in world2['objects']['countries']['geometries'] ]
  
print(country_names[0:5])

['Fiji', 'Tanzania', 'W. Sahara', 'Canada', 'United States of America']

Now let’s check if the names match what is in the Gapminder data set.

import pandas as pd
gapminder = pd.read_json(data.gapminder.url)

common_names = list(set(country_names) & set(gapminder['country']))
missing_names = list(set(gapminder['country']) - set(country_names))
extra_names = list(set(country_names) - set(gapminder['country']))

# names in gapminder and in map data
print("in common:", len(common_names), common_names)

in common: 57 ['Cuba', 'Israel', 'Finland', 'Greece', 'Bangladesh', 'Belgium', 'India', 'Croatia', 'Saudi Arabia', 'Germany', 'Philippines', 'South Africa', 'El Salvador', 'Ecuador', 'Mexico', 'Venezuela', 'Afghanistan', 'United Kingdom', 'Peru', 'Brazil', 'Indonesia', 'Turkey', 'Japan', 'Portugal', 'Rwanda', 'Georgia', 'Egypt', 'Jamaica', 'Pakistan', 'Spain', 'Canada', 'Iran', 'Costa Rica', 'Italy', 'Switzerland', 'Bolivia', 'Iceland', 'Bahamas', 'Lebanon', 'Netherlands', 'Argentina', 'Austria', 'Iraq', 'New Zealand', 'France', 'Colombia', 'Haiti', 'Australia', 'Kenya', 'Chile', 'Ireland', 'South Korea', 'North Korea', 'Norway', 'Poland', 'Nigeria', 'China']

# names in gapminder but not in map data
print("missing in map:", len(missing_names), missing_names)

missing in map: 6 ['Dominican Republic', 'Grenada', 'United States', 'Aruba', 'Hong Kong', 'Barbados']

# names in the map data but not in gapminder
print("extra in map", len(extra_names), extra_names)

extra in map 120 ['Qatar', 'Trinidad and Tobago', 'Algeria', 'Lesotho', 'Gambia', 'Gabon', 'Honduras', 'Nepal', 'Dominican Rep.', 'Cambodia', 'Belize', 'Timor-Leste', 'Sudan', 'Romania', 'Falkland Is.', 'Albania', 'Latvia', 'Tajikistan', 'Estonia', 'Vietnam', 'Burkina Faso', 'Uzbekistan', 'Libya', 'Botswana', 'Russia', 'Tanzania', 'Papua New Guinea', 'S. Sudan', 'Angola', 'Armenia', 'United Arab Emirates', 'Thailand', 'Slovakia', 'Greenland', 'Ghana', 'Guinea', 'Sierra Leone', 'Senegal', 'Eritrea', 'Bhutan', 'Togo', 'Macedonia', 'Belarus', 'Yemen', 'Slovenia', 'Madagascar', 'Vanuatu', 'Mauritania', 'Fr. S. Antarctic Lands', 'Sweden', 'United States of America', 'Dem. Rep. Congo', 'Laos', 'Guyana', 'Bulgaria', 'Nicaragua', 'Uganda', 'Serbia', 'Kazakhstan', 'Paraguay', 'Oman', 'Cyprus', 'Sri Lanka', 'Azerbaijan', 'N. Cyprus', 'Ethiopia', 'Solomon Is.', 'Uruguay', 'Malawi', 'Burundi', 'Bosnia and Herz.', 'Mongolia', 'Kosovo', 'Lithuania', 'Somalia', 'Ukraine', 'Namibia', 'Cameroon', 'Antarctica', 'Tunisia', 'Eq. Guinea', 'Mozambique', 'Czechia', 'Suriname', 'Puerto Rico', 'Hungary', 'Kuwait', 'Zambia', 'Congo', 'Guatemala', 'Niger', 'Somaliland', 'Zimbabwe', 'Central African Rep.', 'Denmark', 'Benin', 'Morocco', 'W. Sahara', 'Chad', 'Guinea-Bissau', 'New Caledonia', 'Jordan', 'Syria', 'Turkmenistan', 'Brunei', 'Mali', 'Kyrgyzstan', 'Taiwan', 'Moldova', 'Panama', 'Malaysia', 'Palestine', 'Fiji', 'Montenegro', "Côte d'Ivoire", 'eSwatini', 'Djibouti', 'Liberia', 'Myanmar', 'Luxembourg']

Chorpoleth of the world

Here is a silly choropleth map of the world that demonstrates how to access the country name.

Note

Tooltip geometry seems to interact badly with these slides. If you open the map in the Vega-Lite editor, the tooltips will appear in the correct place.

vegabrite

Code

world_map <-
  vl_chart() |>
  vl_add_data(url = world2_url, 
               format = list(type = "topojson", feature = "countries")) |>
  vl_calculate("datum.properties.name.length", as = "letters") |>
  vl_calculate("datum.properties.name", as = "country") |>
  vl_mark_geoshape() |>
  vl_encode_fill("letters:Q") |>
  vl_encode_tooltip_array(list("id:N", "country:N", "letters:Q")) |>
  vl_add_properties(width = 800, height = 300)  
world_map

Python/Altair

Code

world = alt.topo_feature(world2_url, feature='countries')
world_map = alt.Chart(world).mark_geoshape().transform_calculate(
  letters = "datum.properties.name.length",
  country = "datum.properties.name"
  ).encode(
    fill = "letters:Q",
    tooltip = ["id:N", "country:N", "letters:Q"]
    )
    
world_map.properties(width = 800, height = 300)

Projection controls

Most projections have additional controls that can be set. For example, some let you choose where to center you map and how much to zoom in.

backwards
upside down
no place like home
interactive

world_map |>
  vl_add_properties(
    width = 800, height = 500,
    projection = list(type = "naturalEarth1", rotate = c(180, 0)))

world_map |>
  vl_add_properties(
    width = 800, height = 500,
    projection = list(type = "naturalEarth1", rotate = c(0, 180)))

world_map |>
  vl_add_properties(
    width = 800, height = 500,
    projection = list(type = "naturalEarth1", rotate = c(85.6681, -42.9634), scale=170))

Code

world_map |>
  vl_add_parameter("rot1", value = 0) |>
  vl_add_parameter("rot2", value = 0) |>
  vl_add_parameter("rot3", value = 0) |>
  vl_add_parameter("scale", value = 200) |>
  vl_bind_range_input("rot1", min = -180, max = 180, step = 1) |>
  vl_bind_range_input("rot2", min = -180, max = 180, step = 1) |>
  vl_bind_range_input("rot3", min = -180, max = 180, step = 1) |>
  vl_bind_range_input("scale", min = 100, max = 500, step = 5) |>
  vl_add_properties(
    width = 800, height = 500,
    projection = list(
      type = "naturalEarth1", 
      scale = list(expr = "scale"),
      rotate = list(expr = "[rot1, rot2, rot3]")
      )
  )

Some useful data

World Factbook Geopolitical Entities, Names, and Codes (GENC)
opendatasoft country codes in various file formats, including GeoJSON
Another source for country names and various standard abbreviations
US states, counties and congressional districts in various file formats, including GeoJSON
Michigan GIS Open Data

More examples

https://nextjournal.com/sdanisch/cartographic-visualization has

Lots of examples
Python/Altair code
Some discussion of projections and how to use them

xkcd and maps

Your turn

Exercise 1 Create a choropleth map of the World using one of the variables in countries, gapminder, or gapminder-health-income from https://cdn.jsdelivr.net/npm/vega-datasets@1.29.0/data/. (Or use some other data of your choosing – but note that you will need country names to match or you will have some extra work to do.)

Exercise 2 Create a choropleth map of the United States where color shows the number of airports in each state.

What other ways might you visualize this information?