Visualizing Spatial Information#

Function geom_livemap() enables a researcher to visualize geospatial information on interactive map.

When building interactive geospatial visualizations with Lets-Plot the visualisation workflow remains the same as it is when building a regular ggplot2 plot.

However, geom_livemap() creates an interactive base-map super-layer and certain limitations do apply comparing to a regular ggplot2 geom-layer:

geom_livemap() must be added as a 1-st layer in plot;
Maximum one geom_livemap() layer is alloed per plot;
Not any type of geometry can be combined with interactive map layer in one plot;
Internet connection to map tiles provider is required.

The following ggplot2 geometry can be used with interactive maps:

geom_point
geom_rect
geom_path
geom_polygon
geom_segment
geom_text
geom_tile
geom_vline, geon_hline
geom_bin2d
geom_contour, geom_contourf
geom_density2d, geom_density2df

from lets_plot import *

LetsPlot.setup_html()

The First Map#

p = ggplot() + ggsize(700, 400)

# Add an empty base-map layer.
p + geom_livemap()

Initial location and zoom level#

The default initial location and zoom level shown on the figure above can be adjusted by using the location and zoom parameters.

# In Lets-Plot coordinates are always encoded as: `longitude` - first and `latitude` - second
p + geom_livemap(location=[29.737809, 60.003151], zoom=11)

Adding ‘geometries’#

`geom_hline()`, `geom_vline()`#

Map’s initial location and zoom parameters are not used very often because map automatically chooses its initial viewport so that all the added geometries are visible.

But the horizontal and vertical lines are exception because they do not have determined bounds. In this case its often necessary to set the desired location and zoom value manually.

lons = [-73 + x for x in range(4)]
lats = [41.5 + 0.5 * y for y in range(4)]

p + geom_livemap(location=[-71.94375, 42.572511], zoom=7)\
  + geom_hline(aes(yintercept=lats), color='red', linetype=2)\
  + geom_vline(aes(xintercept=lons), color='dark_green', linetype=4)

`geom_path()`#

path = {
    'lon': [11.620273, -71.075466, -121.898220, -96.761501], 
    'lat': [48.169211, 42.383651, 37.362240, 32.797371]
}
p + geom_livemap()\
  + geom_path(data=path, mapping=aes(x='lon', y='lat'), size=1, color='magenta')

p + geom_livemap()\
  + geom_path(data=path, mapping=aes(x='lon', y='lat'), size=1, color='magenta')

`geom_point()`#

p + geom_livemap()\
  + geom_path(data=path, mapping=aes(x='lon', y='lat'), size=1, color='magenta')\
  + geom_point(data=path, mapping=aes(x='lon', y='lat'), shape=21, size=7, fill='yellow')

`geom_segment()`#

data = dict(
    x0 = [-97.8,  22.6], y0 = [41.2, -1.3],
    x1 = [ 19.3, -64.5], y1 = [18.7, -15.9],
    Direction = ['east', 'west']
)
p + geom_livemap()\
  + geom_segment(data=data, mapping=aes(x='x0', y='y0', xend='x1', yend='y1', color='Direction'), size=1, linetype=2)

`geom_text()`#

data = {
    'label': ["==    0 ==>", "==   60 ==>", "==  120 ==>", "==  180 ==>", "== -60 ==>", "== -120 ==>"],
    'angle': [0, 60, 120, 180, -60, -120],
}

p + geom_livemap() + geom_point(x=0, y=0, size=140, color='white', alpha=0.6)\
  + geom_text(aes(label = 'label', angle='angle'), data=data,
              x=0, y=0, 
              hjust = "left", vjust = "center",
              size=10, family='monospace', color='red')

`geom_polygon()`#

UK = {
    'lon': [-2.598046, -2.378320, -2.466210, -3.169335, -1.938867, 
            -0.576562, -0.312890, 0.873632, 0.082617, -2.598046], 
    'lat': [51.030349, 51.797754, 53.945750, 54.561879, 55.193929, 
            53.816229, 52.924809, 52.525588, 51.113188, 51.030349]
}
Germany = {
    'lon': [7.685156, 9.926367, 13.661718, 14.101171, 11.464453, 
            12.870703, 8.564062, 8.651953, 6.806250, 6.938085, 7.685156], 
    'lat': [53.294124, 54.049078, 53.608160, 51.305902, 50.221916, 
            48.679365, 48.007575, 49.485266, 50.024691, 51.552493, 53.294124]
}
France = {
    'lon': [-2.246484, 2.367773, 7.245703, 5.268164, 6.586523, 
            2.895117, 2.279882, -0.532617, -0.356835, -2.246484], 
    'lat': [48.095702, 50.586036, 48.795295, 46.365136, 44.169607, 
            43.663114, 43.088157, 43.631315, 46.516550, 48.095702]
}
polygons = dict(
    lon = UK['lon'] + Germany['lon'] + France['lon'],
    lat = UK['lat'] + Germany['lat'] + France['lat'],
    country = ['UK' for _ in UK['lon']] + ['Germany' for _ in Germany['lon']] + ['France' for _ in France['lon']]
)

p + geom_livemap()\
  + geom_polygon(aes('lon', 'lat'),
                 data=polygons, 
                 fill='green', 
                 alpha=0.3)

`geom_rect()`#

bounds = dict(
    lon0 = [min(UK['lon']), min(Germany['lon']), min(France['lon'])],
    lon1 = [max(UK['lon']), max(Germany['lon']), max(France['lon'])],
    lat0 = [min(UK['lat']), min(Germany['lat']), min(France['lat'])],
    lat1 = [max(UK['lat']), max(Germany['lat']), max(France['lat'])],
    country = ['UK', 'Germany', 'France']
)

p + geom_livemap()\
  + geom_rect(aes(xmin='lon0', ymin='lat0', xmax='lon1', ymax='lat1', fill='country'),
              data=bounds, 
              alpha = 0.3) 

Data on Map#

In the examples above we’ve seen how values in input data series can be mapped to aesthetic values on map. For the instance, values in the country series were mapped to the fill color in “geom_rect()” example.

Lets plot another data - the countries population:

pop_UK = 66.65e6
pop_Germany = 83.02e6
pop_France = 66.99e6

Proportional symbols map#

# Create data frame containing countries name, centroid and population.
centroids = dict(
    lon = [(max(xs) - min(xs)) / 2 + min(xs) for xs in [UK['lon'], Germany['lon'], France['lon']]],
    lat = [(max(ys) - min(ys)) / 2 + min(ys)for ys in [UK['lat'], Germany['lat'], France['lat']]],
    country = ['UK', 'Germany', 'France'],
    population = [pop_UK, pop_Germany, pop_France]    
)

p + geom_livemap() + scale_size(range=[10, 50], guide='none')\
  + geom_point(aes('lon', 'lat', color='country', size='population'),
               data=centroids, 
               alpha = 0.8)

Choropleth map#

# Create data frame containing countries name, boundary polygon and population.
polygons_and_pop = polygons.copy()
polygons_and_pop['population'] = [pop_UK for _ in UK['lon']]\
                                + [pop_Germany for _ in Germany['lon']]\
                                + [pop_France for _ in France['lon']]

# Note that in this case it's necessary to define `group` aesthetic.
# Otherwise `geom_polygon` will not be able to split the data into tree "data points".
p + geom_livemap()\
  + geom_polygon(aes('lon', 'lat', fill='population', group='country'), 
                 data=polygons_and_pop, 
                 alpha=0.3)\
  + scale_fill_continuous(format="~s")

Choropleth map (geopandas and `map_join`)#

The dataframe format used in the previouse Choropleth example is not the most convenient way to represent data with geospatial boundaries because each real data-point is duplicated for each vertex in the corresponding polygon.

More convenient would be to store the data-points and their geospatial boundaries in a separate data structures.

In this example let’s use geopandas.GeoDataFrame to store the boundaries.

# Data-points.
population = dict(
    country = ["UK", "Germany", "France"],
    population = [pop_UK, pop_Germany, pop_France]
)

from geopandas import GeoDataFrame, points_from_xy
from shapely.geometry import Polygon, LinearRing

# GeoDataFrame with boundaries.
polygons_gdf = GeoDataFrame(
    data=dict(
        name = ['UK', 'Germany','France']
    ),
    geometry=[
        Polygon(tuple(zip(UK['lon'], UK['lat']))),
        Polygon(tuple(zip(Germany['lon'], Germany['lat']))),
        Polygon(tuple(zip(France['lon'], France['lat'])))
    ]
)

# Use `map` parameter to pass `boundaries` and 
# `map_join` parameter to tell `geom_polygon` how to merge "data" and "map".
p + geom_livemap()\
  + geom_polygon(aes(fill='population'), 
                 data=population, 
                 map=polygons_gdf, 
                 map_join=[['country'],['name']],
                 alpha=0.3)\
  + scale_fill_continuous(format="~s")

Geometries with built-in statistical transformation#

The way these geometries are added to an interactive map is the same as it is done on regular plots.

`geom_bin2d()`#

import numpy as np
np.random.seed(32)
x, y = np.random.multivariate_normal(mean=[0,0], cov=[[1,0],[0,1]], size=50).T
p + geom_livemap() + scale_fill_gradient('white', 'darkgreen')\
  + geom_bin2d(aes(x * 5, y * 5, fill='..density..'), bins=[7, 7], alpha=0.5)

`geom_density2df()`#

p + geom_livemap() + scale_fill_gradient('white', 'darkgreen')\
  + geom_density2df(aes(x * 5, y * 5, fill='..level..'), alpha=0.5)

Visualizing Spatial Information#

The First Map#

Initial location and zoom level#

Adding ‘geometries’#

geom_hline(), geom_vline()#

geom_path()#

geom_point()#

geom_segment()#

geom_text()#

geom_polygon()#

geom_rect()#