Run in Google Colab

View on GitHub

Visualizing Maps¶

Let's draw some maps. 🗺🧐

A dotmap with Altair¶

Let's start with altair.

When your dataset is large, it is nice to enable something called "json data transformer" in altair. What it does is, instead of generating and holding the whole dataset in the memory, to transform the dataset and save into a temporary file. This makes the whole plotting process much more efficient. For more information, check out: https://altair-viz.github.io/user_guide/data_transformers.html

In [1]:

import altair as alt

# saving data into a file rather than embedding into the chart
alt.data_transformers.enable('json')

import altair as alt # saving data into a file rather than embedding into the chart alt.data_transformers.enable('json')

Out[1]:

DataTransformerRegistry.enable('json')

We need a dataset with geographical coordinates. This zipcodes dataset contains the location and zipcode of each zip code area.

In [3]:

from vega_datasets import data

zipcodes_url = data.zipcodes.url
zipcodes = data.zipcodes()
zipcodes.head()

from vega_datasets import data zipcodes_url = data.zipcodes.url zipcodes = data.zipcodes() zipcodes.head()

Out[3]:

	zip_code	latitude	longitude	city	state	county
0	00501	40.922326	-72.637078	Holtsville	NY	Suffolk
1	00544	40.922326	-72.637078	Holtsville	NY	Suffolk
2	00601	18.165273	-66.722583	Adjuntas	PR	Adjuntas
3	00602	18.393103	-67.180953	Aguada	PR	Aguada
4	00603	18.455913	-67.145780	Aguadilla	PR	Aguadilla

In [4]:

zipcodes = data.zipcodes(dtype={'zip_code': 'category'})
zipcodes.head()

zipcodes = data.zipcodes(dtype={'zip_code': 'category'}) zipcodes.head()

Out[4]:

	zip_code	latitude	longitude	city	state	county
0	00501	40.922326	-72.637078	Holtsville	NY	Suffolk
1	00544	40.922326	-72.637078	Holtsville	NY	Suffolk
2	00601	18.165273	-66.722583	Adjuntas	PR	Adjuntas
3	00602	18.393103	-67.180953	Aguada	PR	Aguada
4	00603	18.455913	-67.145780	Aguadilla	PR	Aguadilla

In [5]:

zipcodes.zip_code.dtype

zipcodes.zip_code.dtype

Out[5]:

CategoricalDtype(categories=['00501', '00544', '00601', '00602', '00603', '00604',
                  '00605', '00606', '00610', '00611',
                  ...
                  '99919', '99921', '99922', '99923', '99925', '99926',
                  '99927', '99928', '99929', '99950'],
, ordered=False, categories_dtype=object)

Btw, you'll have fewer issues if you pass URL instead of a dataframe to alt.Chart.

Let's draw it¶

Now we have the dataset loaded and start drawing some plots. Let's say you don't know anything about map projections. What would you try with geographical data? Probably the simplest way is considering (longitude, latitude) as a Cartesian coordinate and directly plot them.

In [6]:

alt.Chart(zipcodes_url).mark_circle().encode(
    x='longitude:Q',
    y='latitude:Q',
)

alt.Chart(zipcodes_url).mark_circle().encode( x='longitude:Q', y='latitude:Q', )

Out[6]:

Actually this itself is a map projection called Equirectangular projection. This projection (or almost a non-projection) is super straight-forward and doesn't require any processing of the data. So, often it is used to just quickly explore geographical data. As you dig deeper, you still want to think about which map projection fits your need best. Don't just use equirectangular projection without any thoughts!

Anyway, let's make it look slighly better by reducing the size of the circles and adjusting the aspect ratio. Q: Can you adjust the width and height?

In [ ]:

alt.Chart(zipcodes_url).mark_circle(
    size=0.5,
).encode(
    x='longitude:Q',
    y='latitude:Q',
    tooltip=['zip_code:N', 'city:N', 'state:N']
).properties(
    width=800,
    height=300
).interactive()

alt.Chart(zipcodes_url).mark_circle( size=0.5, ).encode( x='longitude:Q', y='latitude:Q', tooltip=['zip_code:N', 'city:N', 'state:N'] ).properties( width=800, height=300 ).interactive()

Out[ ]:

But, a much better way to do this is explicitly specifying that they are lat, lng coordinates by using longitude= and latitude=, rather than x= and y=. If you do that, altair automatically adjust the aspect ratio. Q: Can you try it?

In [57]:

alt.Chart(zipcodes_url).mark_circle(
    size=0.5,
).encode(
    longitude='longitude:Q',
    latitude='latitude:Q',
    tooltip=['zip_code:N', 'city:N', 'state:N']
).properties(
    width=800,
    height=300
)

alt.Chart(zipcodes_url).mark_circle( size=0.5, ).encode( longitude='longitude:Q', latitude='latitude:Q', tooltip=['zip_code:N', 'city:N', 'state:N'] ).properties( width=800, height=300 )

Out[57]:

Because the American empire is far-reaching and complicated, the information density of this map is very low (although interesting). Moreover, the US looks twisted because a default projection that is not focused on the US is used.

A common projection for visualizing US data is AlbersUSA, which uses Albers (equal-area) projection. This is a standard projection used in United States Geological Survey and the United States Census Bureau. Albers USA contains a composition of US main land, Alaska, and Hawaii.

To use it, we call project method and specify which variables are longitude and latitude.

Q: use the project method to draw the map in the AlbersUsa projection.

In [59]:

alt.Chart(zipcodes_url).mark_circle(
    size=2,
).encode(
    longitude='longitude:Q',
    latitude='latitude:Q',
    tooltip=['zip_code:N', 'city:N', 'state:N']
).properties(
    width=800,
    height=300
).project(
    type='albersUsa'
    # Supported projections:
    # 'albers', 'albersUsa', 'azimuthalEqualArea', 'azimuthalEquidistant',
    # 'conicConformal', 'conicEqualArea', 'conicEquidistant', 'equalEarth', 'equirectangular', 'gnomonic',
    # 'identity', 'mercator', 'naturalEarth1', 'orthographic', 'stereographic', 'transverseMercator'
).interactive()

alt.Chart(zipcodes_url).mark_circle( size=2, ).encode( longitude='longitude:Q', latitude='latitude:Q', tooltip=['zip_code:N', 'city:N', 'state:N'] ).properties( width=800, height=300 ).project( type='albersUsa' # Supported projections: # 'albers', 'albersUsa', 'azimuthalEqualArea', 'azimuthalEquidistant', # 'conicConformal', 'conicEqualArea', 'conicEquidistant', 'equalEarth', 'equirectangular', 'gnomonic', # 'identity', 'mercator', 'naturalEarth1', 'orthographic', 'stereographic', 'transverseMercator' ).interactive()

Out[59]:

Now we're talking. 😎

Let's visualize the large-scale zipcode patterns. We can use the fact that the zipcodes are hierarchically organized. That is, the first digit captures the largest area divisions and the other digits are about smaller geographical divisions.

Altair provides some data transformation functionalities. One of them is extracting a substring from a variable.

In [22]:

from altair import datum, expr

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate(
    'first_digit', expr.substring(datum.zip_code, 0, 1)
).encode(
    longitude='longitude:Q',
    latitude='latitude:Q',
    color='first_digit:N',
).project(
    type='albersUsa'
).properties(
    width=700,
    height=400,
)

from altair import datum, expr alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate( 'first_digit', expr.substring(datum.zip_code, 0, 1) ).encode( longitude='longitude:Q', latitude='latitude:Q', color='first_digit:N', ).project( type='albersUsa' ).properties( width=700, height=400, )

Out[22]:

For each row (datum), you obtain the zip_code variable and get the substring (imagine Python list slicing), and then you call the result first_digit. Now, you can use this first_digit variable to color the circles. Also note that we specify first_digit as a nominal variable, not quantitative, to obtain a categorical colormap. But we can also play with it too.

Q: Why don't you extract the first two digits, name it as two_digits, and declare that as a quantitative variable? Any interesting patterns? What does it tell us about the history of US?

In [60]:

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate(
    'first_two_digits', expr.substring(datum.zip_code, 0, 2)
).encode(
    longitude='longitude:Q',
    latitude='latitude:Q',
    color=alt.Color('first_two_digits:Q', scale=alt.Scale(scheme='viridis')),
).project(
    type='albersUsa'
).properties(
    width=700,
    height=400,
)

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate( 'first_two_digits', expr.substring(datum.zip_code, 0, 2) ).encode( longitude='longitude:Q', latitude='latitude:Q', color=alt.Color('first_two_digits:Q', scale=alt.Scale(scheme='viridis')), ).project( type='albersUsa' ).properties( width=700, height=400, )

Out[60]:

Q: also try it with declaring the first two digits as a categorical variable

In [11]:

zipcodes

zipcodes

Out[11]:

	zip_code	latitude	longitude	city	state	county
0	00501	40.922326	-72.637078	Holtsville	NY	Suffolk
1	00544	40.922326	-72.637078	Holtsville	NY	Suffolk
2	00601	18.165273	-66.722583	Adjuntas	PR	Adjuntas
3	00602	18.393103	-67.180953	Aguada	PR	Aguada
4	00603	18.455913	-67.145780	Aguadilla	PR	Aguadilla
...	...	...	...	...	...	...
42044	99926	55.094325	-131.566827	Metlakatla	AK	Prince Wales Ketchikan
42045	99927	55.517921	-132.003244	Point Baker	AK	Prince Wales Ketchikan
42046	99928	55.395359	-131.675370	Ward Cove	AK	Ketchikan Gateway
42047	99929	56.449893	-132.364407	Wrangell	AK	Wrangell Petersburg
42048	99950	55.542007	-131.432682	Ketchikan	AK	Ketchikan Gateway

42049 rows × 6 columns

In [ ]:

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate(
    'first_two_digits', expr.substring(datum.zip_code, 0, 2)
).encode(
    longitude='longitude:Q',
    latitude='latitude:Q',
    color='first_two_digits:N',  # note the change here
).project(
    type='albersUsa'
).properties(
    width=700,
    height=400,
)

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate( 'first_two_digits', expr.substring(datum.zip_code, 0, 2) ).encode( longitude='longitude:Q', latitude='latitude:Q', color='first_two_digits:N', # note the change here ).project( type='albersUsa' ).properties( width=700, height=400, )

Out[ ]:

Btw, you can always click "view source" or "open in Vega Editor" to look at the json object that defines this visualization. You can embed this json object on your webpage and easily put up an interactive visualization.

Q: Can you put a tooltip that displays the zipcode when you mouse-over?

In [ ]:

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate(
    'first_two_digits', expr.substring(datum.zip_code, 0, 2)
).encode(
    longitude='longitude:Q',
    latitude='latitude:Q',
    color='first_two_digits:N',
    tooltip=['zip_code:N', 'city:N', 'state:N']  # added tooltip
).project(
    type='albersUsa'
).properties(
    width=700,
    height=400,
)

alt.Chart(zipcodes_url).mark_circle(size=2).transform_calculate( 'first_two_digits', expr.substring(datum.zip_code, 0, 2) ).encode( longitude='longitude:Q', latitude='latitude:Q', color='first_two_digits:N', tooltip=['zip_code:N', 'city:N', 'state:N'] # added tooltip ).project( type='albersUsa' ).properties( width=700, height=400, )

Out[ ]:

In [14]:

zipcodes

zipcodes

Out[14]:

	zip_code	latitude	longitude	city	state	county
0	00501	40.922326	-72.637078	Holtsville	NY	Suffolk
1	00544	40.922326	-72.637078	Holtsville	NY	Suffolk
2	00601	18.165273	-66.722583	Adjuntas	PR	Adjuntas
3	00602	18.393103	-67.180953	Aguada	PR	Aguada
4	00603	18.455913	-67.145780	Aguadilla	PR	Aguadilla
...	...	...	...	...	...	...
42044	99926	55.094325	-131.566827	Metlakatla	AK	Prince Wales Ketchikan
42045	99927	55.517921	-132.003244	Point Baker	AK	Prince Wales Ketchikan
42046	99928	55.395359	-131.675370	Ward Cove	AK	Ketchikan Gateway
42047	99929	56.449893	-132.364407	Wrangell	AK	Wrangell Petersburg
42048	99950	55.542007	-131.432682	Ketchikan	AK	Ketchikan Gateway

42049 rows × 6 columns

Choropleth¶

Let's try some choropleth now. Vega datasets have US county / state boundary data (us_10m) and world country boundary data (world-110m). You can take a look at the boundaries on GitHub (they renders topoJSON files):

• https://github.com/vega/vega-datasets/blob/main/data/us-10m.json
• https://github.com/vega/vega-datasets/blob/main/data/world-110m.json

If you click "Raw" then you can take a look at the actual file, which is hard to read.

Essentially, each file is a large dictionary with the following keys.

In [26]:

usmap = data.us_10m()
usmap.keys()

usmap = data.us_10m() usmap.keys()

Out[26]:

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

In [27]:

usmap['type']

usmap['type']

Out[27]:

'Topology'

In [28]:

usmap['transform']

usmap['transform']

Out[28]:

{'scale': [0.003589294092944858, 0.0005371535195261037],
 'translate': [-179.1473400003406, 17.67439566600018]}

This transformation is used to quantize the data and store the coordinates in integer (easier to store than float type numbers).

https://github.com/topojson/topojson-specification#212-transforms

In [29]:

usmap['objects'].keys()

usmap['objects'].keys()

Out[29]:

dict_keys(['counties', 'states', 'land'])

This data contains not only county-level boundaries (objects) but also states and land boundaries.

In [30]:

usmap['objects']['land']['type'], usmap['objects']['states']['type'], usmap['objects']['counties']['type']

usmap['objects']['land']['type'], usmap['objects']['states']['type'], usmap['objects']['counties']['type']

Out[30]:

('MultiPolygon', 'GeometryCollection', 'GeometryCollection')

land is a multipolygon (one object) and states and counties contains many geometrics (multipolygons) because there are many states (counties). We can look at a state as a set of arcs that define it. It's id captures the identity of the state and is the key to link to other datasets.

In [31]:

state1 = usmap['objects']['states']['geometries'][1]
state1

state1 = usmap['objects']['states']['geometries'][1] state1

Out[31]:

{'type': 'MultiPolygon',
 'arcs': [[[10337]],
  [[10342]],
  [[10341]],
  [[10343]],
  [[10834, 10340]],
  [[10344]],
  [[10345]],
  [[10338]]],
 'id': 15}

The arcs referred here is defined in usmap['arcs'].

In [32]:

usmap['arcs'][:10]

usmap['arcs'][:10]

Out[32]:

[[[15739, 57220], [0, 0]],
 [[15739, 57220], [29, 62], [47, -273]],
 [[15815, 57009], [-6, -86]],
 [[15809, 56923], [0, 0]],
 [[15809, 56923], [-36, -8], [6, -210], [32, 178]],
 [[15811, 56883], [9, -194], [44, -176], [-29, -151], [-24, -319]],
 [[15811, 56043], [-12, -216], [26, -171]],
 [[15825, 55656], [-2, 1]],
 [[15823, 55657], [-19, 10], [26, -424], [-26, -52]],
 [[15804, 55191], [-30, -72], [-47, -344]]]

It seems pretty daunting to work with this dataset, right? But fortunately people have already built tools to handle such data.

In [33]:

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

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

In [34]:

states

states

Out[34]:

UrlData({
  format: TopoDataFormat({
    feature: 'states',
    type: 'topojson'
  }),
  url: 'https://cdn.jsdelivr.net/npm/vega-datasets@v1.29.0/data/us-10m.json'
})

Can you find a mark for geographical shapes from here https://altair-viz.github.io/user_guide/marks/index.html# and draw the states?

In [37]:

# Can you find a mark for geographical shapes from here https://altair-viz.github.io/user_guide/marks/index.html# and draw the states?
alt.Chart(states).mark_geoshape(
    fill='lightgray',
    stroke='black'
).encode(
    tooltip='id:N'
).properties(
    width=800,
    height=400
)

# Can you find a mark for geographical shapes from here https://altair-viz.github.io/user_guide/marks/index.html# and draw the states? alt.Chart(states).mark_geoshape( fill='lightgray', stroke='black' ).encode( tooltip='id:N' ).properties( width=800, height=400 )

Out[37]:

And then project it using the albersUsa?

In [39]:

alt.Chart(states).mark_geoshape(
    fill='lightgray',
    stroke='black'
).encode(
    tooltip='id:N'
).properties(
    width=800,
    height=400
).project(
    type='albersUsa'
)

alt.Chart(states).mark_geoshape( fill='lightgray', stroke='black' ).encode( tooltip='id:N' ).properties( width=800, height=400 ).project( type='albersUsa' )

Out[39]:

Can you do the same thing with counties and draw county boundaries?

In [62]:

counties = alt.topo_feature(data.us_10m.url, 'counties')
alt.Chart(counties).mark_geoshape(
    fill='lightgray',
    stroke='black'
).encode(
    tooltip='id:N'
).properties(
    width=800,
    height=400
).project(
    type='albersUsa'
)

counties = alt.topo_feature(data.us_10m.url, 'counties') alt.Chart(counties).mark_geoshape( fill='lightgray', stroke='black' ).encode( tooltip='id:N' ).properties( width=800, height=400 ).project( type='albersUsa' )

Out[62]:

Let's load some county-level unemployment data.

In [41]:

unemp_data = data.unemployment(sep='\t')
unemp_data.head()

unemp_data = data.unemployment(sep='\t') unemp_data.head()

Out[41]:

	id	rate
0	1001	0.097
1	1003	0.091
2	1005	0.134
3	1007	0.121
4	1009	0.099

This dataset has unemployment rate. When? I don't know. We don't care about data provenance here because the goal is quickly trying out choropleth. But if you're working with a real dataset, you should be very sensitive about the provenance of your dataset. Make sure you understand where the data came from and how it was processed.

Anyway, for each county specified with id. To combine two datasets, we use "Lookup transform" - https://vega.github.io/vega/docs/transforms/lookup/. Essentially, we use the id in the map data to look up (again) id field in the unemp_data and then bring in the rate variable. Then, we can use that rate variable to encode the color of the geoshape mark.

In [42]:

alt.Chart(counties).mark_geoshape().project(
    type='albersUsa'
).transform_lookup(
    lookup='id',
    from_=alt.LookupData(unemp_data, 'id', ['rate'])
).encode(
    color='rate:Q'
).properties(
    width=700,
    height=400
)

alt.Chart(counties).mark_geoshape().project( type='albersUsa' ).transform_lookup( lookup='id', from_=alt.LookupData(unemp_data, 'id', ['rate']) ).encode( color='rate:Q' ).properties( width=700, height=400 )

Out[42]:

There you have it, a nice choropleth map. 😎

Leaflet¶

Another useful tool is Leaflet. It allows you to use various map tile data (Google maps, Open streetmap, ...) with many types of marks (points, heatmap, etc.). Leaflet.js is one of the easiest options to do that on the web, and there is a Python bridge of it: https://github.com/jupyter-widgets/ipyleaflet.

You can install it simply by

pip install ipyleaflet

It is quite easy to display an interactive map with it. You can also add markers, polygons, etc.

In [48]:

from ipyleaflet import Map, Marker

center = (39.1737568,-86.5220677)

m = Map(center=center, zoom=15)

marker = Marker(location=center, draggable=False)
m.add(marker)

display(m)

marker.location = center

from ipyleaflet import Map, Marker center = (39.1737568,-86.5220677) m = Map(center=center, zoom=15) marker = Marker(location=center, draggable=False) m.add(marker) display(m) marker.location = center

Map(center=[39.1737568, -86.5220677], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_titl…

ipyleaflet lets you use various basemaps. You can find the list of available basemaps here: https://leaflet-extras.github.io/leaflet-providers/preview/

For instance, let's use OpenStreetMap and zoom into the Luddy School building.

In [49]:

from ipyleaflet import basemaps

center = (39.172681590059604, -86.5233788123735)
zoom = 17

Map(basemap=basemaps.OpenStreetMap.HOT, center=center, zoom=zoom)

from ipyleaflet import basemaps center = (39.172681590059604, -86.5233788123735) zoom = 17 Map(basemap=basemaps.OpenStreetMap.HOT, center=center, zoom=zoom)

Out[49]:

Map(center=[39.172681590059604, -86.5233788123735], controls=(ZoomControl(options=['position', 'zoom_in_text',…

Let's create a heatmap. First create a small list of random points around the center.

In [50]:

center

center

Out[50]:

(39.172681590059604, -86.5233788123735)

In [52]:

from random import uniform

random_points = [
    [
        uniform(center[0]-0.01, center[0]+0.01), 
        uniform(center[1]-0.01, center[1]+0.01), 
        uniform(0, 5)
    ] for x in range(100)
]
random_points[:2]

from random import uniform random_points = [ [ uniform(center[0]-0.01, center[0]+0.01), uniform(center[1]-0.01, center[1]+0.01), uniform(0, 5) ] for x in range(100) ] random_points[:2]

Out[52]:

[[39.1668122403533, -86.52327515084703, 2.1730825811477006],
 [39.16729783526083, -86.52370948019752, 4.741391103458917]]

Q: using these random points, can you create a heatmap?

documentation: https://ipyleaflet.readthedocs.io/en/latest/layers/heatmap.html

In [56]:

from ipyleaflet import Map, Heatmap

m = Map(center=center, zoom=zoom)
heatmap = Heatmap(
    locations=random_points,
    radius=20,
    blur=15,
    max_zoom=1,
    gradient={
        0.0: 'white',
        1.0: 'red'
    }
)
m.add_layer(heatmap)
m

from ipyleaflet import Map, Heatmap m = Map(center=center, zoom=zoom) heatmap = Heatmap( locations=random_points, radius=20, blur=15, max_zoom=1, gradient={ 0.0: 'white', 1.0: 'red' } ) m.add_layer(heatmap) m

Out[56]:

Map(center=[39.172681590059604, -86.5233788123735], controls=(ZoomControl(options=['position', 'zoom_in_text',…