Getting a tile image in your desired dimensions from MapTiles.jl and TileProviders.jl

tiles_as_raster.jl

#=
# Getting (noninteractive) raster images from the MapTiles/TileProviders ecosystem
=#

import HTTP # to run the actual request and fetch the tiles
import ImageMagick # to read returned images, since they may be any bitmap format
using MapTiles, TileProviders # Julia map-tile infrastructure
# Julia geographic infrastructure
using Extents 
import GeoInterface as GI
import GeoFormatTypes as GFT
# import Proj # reprojection
import Rasters # just to have a nice interface somehow


bbox = Extent(X = (-1, 5.65), Y = (-1, 5.65))
provider = TileProviders.OpenStreetMap()
#=
The API is simply
```julia
Rasters.Raster(TileProviders.OpenStreetMap(), Extents.Extent(X = (1, 2), Y = (3, 4)); size = (1000, 1000))`.
```
=#
# Stolen from Tyler.jl
function z_index(extent::Union{Rect,Extent}, res::NamedTuple, crs)
    # Calculate the number of tiles at each z and get the one
    # closest to the resolution `res`
    target_ntiles = prod(map(r -> r / 256, res))
    tiles_at_z = map(1:24) do z
        length(TileGrid(extent, z, crs))
    end
    return findmin(x -> abs(x - target_ntiles), tiles_at_z)[2]
end


function Rasters.Raster(provider::TileProviders.Provider, bbox::Extent; res = nothing, size = nothing, min_zoom_level = 0, max_zoom_level = 16)
    # First, handle keyword arguments
    @assert isnothing(size) || isnothing(res) "You must provide either `size` or `res`, but not both."
    _size = if isnothing(size) 
        # convert resolution to size using bbox and round(Int, x)
        (round(Int, (bbox.X[2] - bbox.X[1]) / res[1]), round(Int, (bbox.Y[2] - bbox.Y[1]) / res[2]))
    else
        (first(size), last(size))
    end
    # Obtain the optimal Z-index that covers the bbox at the desired resolution.
    optimal_z_index = clamp(z_index(bbox, (X=_size[2], Y=_size[1]), MapTiles.WGS84()), min_zoom_level, max_zoom_level)
    # Generate a `TileGrid` from our zoom level and bbox.
    tilegrid = MapTiles.TileGrid(bbox, optimal_z_index, MapTiles.WGS84())
    # Compute the dimensions of the tile grid, so we can feed them into a 
    # Raster later.
    tilegrid_extent = Extents.extent(tilegrid, MapTiles.WGS84())
    tilegrid_size = tile_widths .* length.(tilegrid.grid.indices)
    # We need to know the start and end indices of the tile grid, so we can 
    # place the tiles in the right place.
    tile_start_idxs = minimum(first.(Tuple.(tilegrid.grid))), minimum(last.(Tuple.(tilegrid.grid)))
    tile_end_idxs = maximum(first.(Tuple.(tilegrid.grid))), maximum(last.(Tuple.(tilegrid.grid)))
    # Using the size information, we initiate an `RGBA{Float32}` image array.
    # You can later convert to whichever size / type you want by simply broadcasting.
    image_receptacle = fill(RGBAf(0,0,0,1), tilegrid_size)
    # Now, we iterate over the tiles, and read and then place them into the array.
    for tile in tilegrid
        # Download the tile
        url = TileProviders.geturl(provider, tile.x, tile.y, tile.z)
        result = HTTP.get(url)
        # Read into an in-memory array (Images.jl layout)
        img = ImageMagick.readblob(result.body)
        # The thing with the y indices is that they go in the reverse of the natural order.
        # So, we simply subtract the y index from the end index to get the correct placement.
        image_start_relative = (
            (tile.x) - tile_start_idxs[1], 
            tile_end_idxs[2] - tile.y,
        )
        # The absolute start is simply the relative start times the tile width.
        image_start_absolute = (image_start_relative .* tile_widths)
        # The indices for the view into the receptacle are the absolute start 
        # plus one, to the absolute end.
        idxs = (:).(image_start_absolute .+ 1, image_start_absolute .+ tile_widths)
        @debug image_start_relative image_start_absolute idxs
        # Place the tile into the receptacle.  Note that we rotate the image to 
        # be in the correct orientation.
        image_receptacle[idxs...] .= rotr90(img) # change to Julia memory layout
    end
    # Now, we have a complete image.  We can create a `Raster` from it.
    ras = Raster(
        image_receptacle, 
        ( # dims
            Rasters.X(LinRange(tilegrid_extent.X..., tilegrid_size[1])), 
            Rasters.Y(LinRange(tilegrid_extent.Y..., tilegrid_size[2]))
        );
        crs = MapTiles.WGS84(), # TODO: should this be a different CRS?
    )
    # image(ras; axis = (; aspect = DataAspect()))
    return ras
end

This is now implemented as Tyler.basemap in a PR to Tyler.