Географические фигуры — Geoshape

mark_geoshape представляет произвольные фигуры, геометрия которых определяется заданными пространственными данными.

Свойства географических фигур

Примитив geoshape поддерживает любые стандартные свойства.

Простая карта

Altair работает с разными географическими форматами, включая GeoJSON и TopoJSON. Часто удобнее всего использовать GeoDataFrame. Загрузим набор Natural Earth 110m Cultural Vectors и создадим простую карту:

import altair as alt
from altair.datasets import data
import geopandas as gpd

url = "https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip"
gdf_ne = gpd.read_file(url)  # zipped shapefile
gdf_ne = gdf_ne[["NAME", "CONTINENT", "POP_EST", 'geometry']]

alt.Chart(gdf_ne).mark_geoshape()

Здесь Altair применяет стандартную синюю заливку fill и проекцию equalEarth. Цвета и толщину границ можно настроить стандартными свойствами, а проекцию — явно задать методом project:

alt.Chart(gdf_ne).mark_geoshape(
    fill='lightgrey', stroke='white', strokeWidth=0.5
).project(
    type='albers'
)

Область показа и фильтрация

По умолчанию Altair автоматически подбирает проекцию, чтобы все данные помещались в ширину и высоту диаграммы. Выделить нужный регион можно несколькими способами:

1. Отфильтровать данные внутри GeoDataFrame.
2. Применить transform_filter.
3. Указать scale для масштаба и translate для перемещения в project.
4. Указать область fit в project и clip=True в свойствах примитива.

Применим эти подходы для показа континентальной Африки.

1. Фильтрация в GeoDataFrame:

gdf_sel = gdf_ne.query("CONTINENT == 'Africa'")

alt.Chart(gdf_sel).mark_geoshape()

2. Фильтрация через transform_filter:

alt.Chart(gdf_ne).mark_geoshape().transform_filter(
    alt.datum.CONTINENT == 'Africa'
)

3. Масштаб и перемещение в project:

alt.Chart(gdf_ne).mark_geoshape().project(
    scale=200,
    translate=[160, 160]  # lon, lat
)

4. Область fit в project и clip=True:

extent_roi = gdf_ne.query("CONTINENT == 'Africa'")
xmin, ymin, xmax, ymax = extent_roi.total_bounds

# fit object should be a GeoJSON-like Feature or FeatureCollection 
extent_roi_feature = {
    "type": "Feature", 
    "geometry": {"type": "Polygon", 
                 "coordinates": [[
                     [xmax, ymax],
                     [xmax, ymin],
                     [xmin, ymin],
                     [xmin, ymax],
                     [xmax, ymax]]]},
    "properties": {}
}

alt.Chart(gdf_ne).mark_geoshape(clip=True).project(
    fit=extent_roi_feature
)

Декартовы координаты

Стандартная проекция equalEarth правильно передаёт относительные площади суши. Она предполагает геометрию в градусах с координатами долготы и широты. Другая распространённая система — двумерные декартовы координаты, не учитывающие кривизну Земли.

В следующем примере геометрия отображается непосредственно в исходных координатах через проекцию identity. Также нужен reflectY, поскольку в Canvas и SVG положительная ось y направлена вниз:

alt.Chart(gdf_sel).mark_geoshape().project(
    type='identity',
    reflectY=True
)

Отображение полигонов

Сопоставим столбец NAME с цветом:

alt.Chart(gdf_sel).mark_geoshape().encode(
    color='NAME:N'
)

Каждая страна представлена полигоном или мультиполигоном, поэтому заливку и обводку можно задать отдельно:

alt.Chart(gdf_sel).mark_geoshape(
    stroke='white',
    strokeWidth=1.5
).encode(
    fill='NAME:N'
)

Отображение линий

По умолчанию цвет mark_geoshape используется для заливки, а не обводки. Если источник содержит линии или мультилинии, явно задайте filled=False.

Сравните:

gs_line = gpd.GeoSeries.from_wkt(['LINESTRING (0 0, 1 1, 0 2, 2 2, -1 1, 1 0)'])
alt.Chart(gs_line).mark_geoshape().project(
    type='identity',
    reflectY=True
)

с:

gs_line = gpd.GeoSeries.from_wkt(['LINESTRING (0 0, 1 1, 0 2, 2 2, -1 1, 1 0)'])
alt.Chart(gs_line).mark_geoshape(
    filled=False
).project(
    type='identity',
    reflectY=True
)

Таким способом полигоны также можно оформить как линейные геометрии:

alt.Chart(gdf_sel).mark_geoshape(
    filled=False,
    strokeWidth=1.5
).encode(
    stroke='NAME:N'
)

Отображение точек

mark_geoshape рисует точки, если они заданы геометриями Point в GeoDataFrame. Сначала заменим полигоны их центроидами и покажем результат:

# .copy() to prevent changing the original `gdf_sel` variable
# derive centroid in a projected CRS (in meters) and visualize in a geographic CRS (in degrees).
gdf_centroid = gpd.GeoDataFrame(
    data=gdf_sel.copy(),
    geometry=gdf_sel.geometry.to_crs(epsg=3857).centroid.to_crs(epsg=4326)
)

alt.Chart(gdf_centroid).mark_geoshape()

Для кодирования размера этих точек нужно использовать mark_circle вместе с каналами latitude и longitude:

gdf_centroid["lon"] = gdf_centroid.geometry.x
gdf_centroid["lat"] = gdf_centroid.geometry.y

alt.Chart(gdf_centroid).mark_circle().encode(
    longitude="lon:Q", latitude="lat:Q", size="POP_EST:Q"
)

В Altair есть выражения для географических объектов. Например, центроиды можно вычислить через geoCentroid:

basemap = alt.Chart(gdf_sel).mark_geoshape(
     fill='lightgray', stroke='white', strokeWidth=0.5
)

bubbles = alt.Chart(gdf_sel).transform_calculate(
    centroid=alt.expr.geoCentroid(None, alt.datum)
).mark_circle(
    stroke='black'
).encode(
    longitude='centroid[0]:Q',
    latitude='centroid[1]:Q',
    size="POP_EST:Q"
)

(basemap + bubbles).project(
    type='identity', reflectY=True
)

Картограммы

Вместо пузырьков с численностью населения можно построить картограмму: географическую тепловую карту, где цвет региона определяется значением столбца DataFrame.

alt.Chart(gdf_sel).mark_geoshape().encode(
    color='POP_EST'
)

При чтении картограмм нужна осторожность: цвет отражает нужное значение, но размер фигуры определяется площадью страны. Важные значения небольших стран легко пропустить, потому что их плохо видно, например при отображении плотности населения.

Соединение наборов данных

Иногда данные разделены: обычный DataFrame содержит значения, а GeoDataFrame — геометрию. Преобразование lookup позволяет получить связанную информацию из второго набора.

Возможны два направления:

1. Использовать GeoDataFrame с геометриями как источник и добавить значения другого DataFrame.
2. Использовать DataFrame как источник и добавить геометрии GeoDataFrame.

Выбор зависит от задачи. Сначала найдём поле rate в df_us_unemp, используя gdf_us_counties как источник:

import altair as alt
from altair.datasets import data

gdf_us_counties = data.us_10m(layer="counties")
df_us_unemp = data.unemployment()

alt.Chart(gdf_us_counties).mark_geoshape().transform_lookup(
    lookup='id',
    from_=alt.LookupData(data=df_us_unemp, key='id', fields=['rate'])
).encode(
    alt.Color('rate:Q')
).project(
    type='albersUsa'
)

Затем обратный подход: найдём geometry и type в gdf_us_counties, используя df_us_unemp как источник:

alt.Chart(df_us_unemp).mark_geoshape().transform_lookup(
    lookup='id',
    from_=alt.LookupData(data=gdf_us_counties, key='id', fields=['geometry', 'type'])
).encode(
    alt.Color('rate:Q')
).project(
    type='albersUsa'
)

Классификация для картограмм

Картограммы могут показывать не только непрерывные значения, но и дискретные уровни. Искусственные группы при дискретизации требуют осторожности, однако естественные пороги полезно выделить явно.

Определим вспомогательную функцию classify() для сравнения подходов. Начнём с картограммы безработицы по округам США за 2018 год с линейной шкалой:

import altair as alt
from altair.datasets import data

def classify(type, domain=None, nice=False, title=None):
    # define data
    us_counties = alt.topo_feature(data.us_10m.url, "counties")
    us_unemp = data.unemployment.url

    # define choropleth scale
    if "threshold" in type:
        scale = alt.Scale(type=type, domain=domain, scheme="inferno")
    else:
        scale = alt.Scale(type=type, nice=nice, scheme="inferno")

    # define title
    if title is None:
        title = type

    # define choropleth chart
    choropleth = (
        alt.Chart(us_counties, title=title)
        .mark_geoshape()
        .transform_lookup(
            lookup="id", from_=alt.LookupData(data=us_unemp, key="id", fields=["rate"])
        )
        .encode(
            alt.Color(
                "rate:Q",
                scale=scale,
                legend=alt.Legend(
                    direction="horizontal", orient="bottom", format=".1%"
                ),
            )
        )
        .project(type="albersUsa")
    )
    return choropleth

classify(type='linear')

Уровень безработицы показан в процентах на карте с mark_geoshape(). При линейной шкале каждое значение округа получает свой цвет. Это полезно, но часто распределение удобнее сгруппировать в классы, присваивая всем значениям и геометриям класса один цвет.

Рассмотрим несколько способов.

quantile делит набор (domain) на интервалы примерно одинаковой численности. В каждом классе примерно столько же значений или геометрий — equal counts:

alt.Scale(type='quantile')

Применим в нашей функции:

classify(type='quantile', title=['quantile', 'equal counts'])

quantize делит диапазон данных на равные интервалы. Число значений в классах может различаться, но шаг одинаков — equal range:

alt.Scale(type='quantize')

Применим в функции:

classify(type='quantize', title=['quantize', 'equal range'])

quantize можно сочетать с nice, чтобы перед дискретизацией привести границы диапазона к удобным значениям:

alt.Scale(type='quantize', nice=True)

Применим в функции:

classify(type='quantize', nice=True, title=['quantize', 'equal range nice'])

threshold делит данные на классы с вручную заданными порогами:

alt.Scale(type='threshold', domain=[0.05, 0.20])

Применим в функции:

classify(type='threshold', domain=[0.05, 0.20])

Получатся три класса: ниже 0.05, от 0.05 до 0.20 и от 0.20 выше.

Какой способ классификации оптимален? Как обычно, это зависит от задачи.

Есть популярный метод поиска естественных границ (natural breaks): он повышает сходство значений внутри класса и различия между классами. Это алгоритм Фишера — Дженкса, похожий на одномерный k-means.

Пакет Python jenkspy позволяет вычислить такие границы:

>>> from jenkspy import JenksNaturalBreaks
>>> jnb = JenksNaturalBreaks(5)
>>> jnb.fit(df_us_unemp['rate'])
>>> jnb.inner_breaks_
[0.061, 0.088, 0.116, 0.161]

Применим их в функции:

classify(type='threshold', domain=[0.061, 0.088, 0.116, 0.161],
        title=['threshold Jenks','natural breaks'])

Особенности:

• У шкал quantize и quantile по умолчанию пять классов. Число меняется через SchemeParams(). Например, вместо scheme='turbo' можно задать scheme=alt.SchemeParams('turbo', count=2) для двух классов.
• Метод естественных границ находит оптимальные пороги при заданном числе классов. Для выбора этого числа можно исследовать показатель качества соответствия дисперсии (GVF), также называемый методом оптимизации Дженкса.

Повторение карты

RepeatChart, создаваемый методом Chart.repeat, удобен для горизонтального или вертикального объединения представлений многомерного набора.

В следующем примере используются столбцы population, engineers и hurricanes для штатов США. Геометрия states создаётся функцией topo_feature с параметрами url и feature, извлекающей объекты из TopoJSON по URL.

Столбцы передаются списком в .repeat(), а в цветовом кодировании используется alt.repeat('row'):

import altair as alt
from altair.datasets import data

states = alt.topo_feature(data.us_10m.url, 'states')
source = data.population_engineers_hurricanes.url
variable_list = ['population', 'engineers', 'hurricanes']

alt.Chart(states).mark_geoshape(tooltip=True).encode(
    alt.Color(alt.repeat('row'), type='quantitative')
).transform_lookup(
    lookup='id',
    from_=alt.LookupData(source, 'id', variable_list)
).project(
    type='albersUsa'
).repeat(
    row=variable_list
).resolve_scale(
    color='independent'
)

Разбиение карты на панели

FacetChart, создаваемый методом Chart.facet, удобен для объединения представлений, когда одно поле содержит несколько переменных.

Пока issue #2369 не решён, обычное разбиение не работает для географических визуализаций:

source = data.population_engineers_hurricanes().melt(id_vars=['state', 'id'])
us_states = data.us_10m(layer="states")
gdf_comb = gpd.GeoDataFrame(source.join(us_states, on='id', rsuffix='_y'))

alt.Chart(gdf_comb).mark_geoshape().encode(
    color=alt.Color('value:Q'),
    facet=alt.Facet('variable:N').columns(3)
).properties(
    width=180,
    height=130
).resolve_scale('independent')

Есть два обходных подхода. Можно передать географические данные через transform_lookup, а не через Chart, как в примере групп доходов по штатам. Либо отфильтровать данные вручную в pandas и объединить отдельные карты, как здесь:

source = data.population_engineers_hurricanes().melt(id_vars=['state', 'id'])
us_states = data.us_10m(layer="states")
gdf_comb = gpd.GeoDataFrame(source.join(us_states, on='id', rsuffix='_y'))

alt.concat(
    *(
        alt.Chart(gdf_comb[gdf_comb.variable == var], title=var)
        .mark_geoshape()
        .encode(
            color=alt.Color(
                "value:Q", legend=alt.Legend(orient="bottom", direction="horizontal")
            )
        )
        .project('albersUsa')
        .properties(width=180, height=130)
        for var in gdf_comb.variable.unique()
    ),
    columns=3
).resolve_scale(color="independent")

Интерактивность

Карта часто используется вместе с другой диаграммой. Покажем интерактивную карту штатов США и столбчатую диаграмму 15 самых населённых штатов. alt.selection_point() задаёт выбор по щелчку левой кнопкой:

import altair as alt
from altair.datasets import data

# load the data
us_states = data.us_10m(layer="states")
us_population = data.population_engineers_hurricanes()[["state", "id", "population"]]

# define a pointer selection
click_state = alt.selection_point(fields=["state"])
# define a condition on the opacity encoding depending on the selection
opacity = alt.when(click_state).then(alt.value(1)).otherwise(alt.value(0.2))

# create a choropleth map using a lookup transform
choropleth = (
    alt.Chart(us_states)
    .mark_geoshape()
    .transform_lookup(
        lookup="id", from_=alt.LookupData(us_population, "id", ["population", "state"])
    )
    .encode(
        color="population:Q",
        opacity=opacity,
        tooltip=["state:N", "population:Q"],
    )
    .project(type="albersUsa")
)

# create a bar chart with the same conditional ``opacity`` encoding.
bars = (
    alt.Chart(
        us_population.nlargest(15, "population"), title="Top 15 states by population"
    )
    .mark_bar()
    .encode(
        x="population",
        opacity=opacity,
        color="population",
        y=alt.Y("state").sort("-x"),
    )
)

(choropleth & bars).add_params(click_state)

Связь двусторонняя: щелчок по геометрии или столбцу, в том числе с Shift для множественного выбора, задаёт выбранным элементам opacity=1, а остальным — opacity=0.2.

Возможен и интервальный выбор, как в примере карты землетрясений.

Выражения

Выражения Altair работают и в географической визуализации. Здесь землетрясения показаны на глобусе в проекции orthographic, а Землю можно вращать по одной оси (rotate0). Функция sphere создаёт диск Земли на фоне.

Геометрия точек землетрясений имеет координаты XYZ: долгота, широта и глубина. Координаты извлекаются напрямую из вложенного поля geometry для отображения кругов, без предварительного создания трёх отдельных столбцов:

import altair as alt
from altair.datasets import data

# load data
gdf_quakies = data.earthquakes()
gdf_world = data.world_110m(layer="countries")

# define parameters
range0 = alt.binding_range(min=-180, max=180, step=5, name='rotate longitude ')
rotate0 = alt.param(value=120, bind=range0)
hover = alt.selection_point(on="pointerover", clear="pointerout")

# world disk
sphere = alt.Chart(alt.sphere()).mark_geoshape(
    fill="aliceblue", stroke="black", strokeWidth=1.5
)

# countries as shapes
world = alt.Chart(gdf_world).mark_geoshape(
    fill="mintcream", stroke="black", strokeWidth=0.35
)

# earthquakes as circles with fill for depth and size for magnitude
# the hover param is added on the mar_circle only
quakes = (
    alt.Chart(gdf_quakies)
    .mark_circle(opacity=0.35, tooltip=True, stroke="black")
    .transform_calculate(
        lon="datum.geometry.coordinates[0]",
        lat="datum.geometry.coordinates[1]",
        depth="datum.geometry.coordinates[2]",
    )
    .transform_filter(
        ((rotate0 * -1 - 90 < alt.datum.lon) & (alt.datum.lon < rotate0 * -1 + 90)).expr
    )
    .encode(
        longitude="lon:Q",
        latitude="lat:Q",
        strokeWidth=alt.when(hover, empty=False).then(alt.value(1)).otherwise(alt.value(0)),
        size=alt.Size(
            "mag:Q",
            scale=alt.Scale(type="pow", range=[1, 1000], domain=[0, 6], exponent=4),
        ),
        fill=alt.Fill(
            "depth:Q", scale=alt.Scale(scheme="lightorange", domain=[0, 400])
        ),
    )
    .add_params(hover, rotate0)
)

# define projection and add the rotation param for all layers
comb = alt.layer(sphere, world, quakes).project(
    type="orthographic",
    rotate=alt.expr(f"[{rotate0.name}, 0, 0]")
)
comb

Землетрясения отображаются кругами и фильтруются, когда выходят за видимую часть глобуса. Выделение при наведении позволяет подробнее изучить отдельное событие.

Тайловые карты

Чтобы использовать тайловые карты, например OpenStreetMap, как фон mark_geoshape, можно сочетать Altair с пакетом Altair Tiles.