altair_osm_tiles.py

import altair as alt
from vega_datasets import data

source = alt.topo_feature(data.world_110m.url, "countries")
osm_url, otm_url = ('https://tile.openstreetmap.org/', 'https://tile.opentopomap.org/')
select_urls = alt.binding_select(options=[osm_url, otm_url], name='select tile service')
param_urls = alt.param(bind=select_urls, value=osm_url)
param_tx = alt.param(expr="width / 2")
param_ty = alt.param(expr="height / 2")
param_base_tile_size = alt.param(value=256)
range_z = alt.binding_range(min=2, max=13, step=0.05, name="zoom level")
param_z = alt.param(value=2.75, bind=range_z)
range_x = alt.binding_range(min=-180, max=180, step=0.05, name="rotate longitude")
param_x = alt.param(value=-5.9025, bind=range_x)
range_y = alt.binding_range(min=-60, max=60, step=0.05, name="center latitude")
param_y = alt.param(value=52.56, bind=range_y)
param_tile_url = alt.param(expr=f'{param_urls.name}')
param_zoom = alt.param(expr=f"ceil({param_z.name})")
param_tiles_count = alt.param(expr=f"pow(2, {param_zoom.name})")
param_tile_size = alt.param(
    expr=f"{param_base_tile_size.name} * pow(2, {param_z.name} - {param_zoom.name})"
)
param_base_point = alt.param(expr=f"invert('projection', [0, 0])")
param_dii = alt.param(
    expr=f"({param_base_point.name}[0] + 180) / 360 * {param_tiles_count.name}"
)
param_di = alt.param(expr=f"floor({param_dii.name})")
param_dx = alt.param(
    expr=f"round((floor({param_dii.name}) - {param_dii.name}) * {param_tile_size.name})"
)
param_djj = alt.param(
    expr=f"(1 - log(tan({param_base_point.name}[1] * PI / 180) + 1 / cos({param_base_point.name}[1] * PI / 180)) / PI) / 2 * {param_tiles_count.name}"
)
param_dj = alt.param(expr=f"floor({param_djj.name})")
param_dy = alt.param(
    expr=f"round((floor({param_djj.name})-{param_djj.name}) * {param_tile_size.name})"
)
tile_list = alt.sequence(0, 4, as_="a", name="tile_list")


image_tiles = (
    alt.Chart(tile_list)
    .mark_image(
        width=alt.expr(f"{param_tile_size.name}"),
        height=alt.expr(f"{param_tile_size.name}"),
        clip=True,
    )
    .transform_calculate(b=f"sequence(0, 4)")
    .transform_flatten(["b"])
    .transform_calculate(
        url=f"{param_tile_url.name} + {param_zoom.name} + '/' + (datum.a + {param_di.name} + {param_tiles_count.name}) % {param_tiles_count.name} + '/' + ((datum.b + {param_dj.name})) + '.png'",
        x=f"(datum.a * {param_tile_size.name} + {param_dx.name}) + ({param_tile_size.name} / 2)",
        y=f"(datum.b * {param_tile_size.name} + {param_dy.name}) + ({param_tile_size.name} / 2)",
    )
    .encode(
        x=alt.X("x:Q").scale(None), y=alt.Y("y:Q").scale(None), url=alt.Url("url:N")
    )
)

geoshape_countries = (
    alt.Chart(source, width=400, height=400)
    .mark_geoshape(
        stroke="orange", 
        strokeWidth=2, 
        fillOpacity=0.1
    )
    .encode(fill="id:Q")
    .project(
        type="mercator",
        scale=alt.expr(
            f"{param_base_tile_size.name} * pow(2, {param_z.name}) / (2 * PI)"
        ),
        rotate=alt.expr(f"[{param_x.name}, 0, 0]"),
        center=alt.expr(f"[0, {param_y.name}]"),
        translate=alt.expr(f"[{param_tx.name}, {param_ty.name}]"),
    )
)

text_attrib = alt.Chart().mark_text(
    text='(C) OpenStreetMap contributors',
    dx=-85,
    dy=-10
).encode(
    x=alt.value(alt.expr('width')), 
    y=alt.value(alt.expr('height'))
)

alt.layer(
    image_tiles,
    geoshape_countries,
    text_attrib
).add_params(
    param_urls,
    param_tile_url,
    param_zoom,
    param_tiles_count,
    param_tile_size,
    param_base_point,
    param_dii,
    param_di,
    param_dx,
    param_djj,
    param_dj,
    param_dy,
    param_y,
    param_x,
    param_z,
    param_tx,
    param_ty,
    param_base_tile_size,
)