Add a function basemap to get a static image from Tyler, and some converts for image

docs/Project.toml

@@ -10,11 +10,14 @@ DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
+Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LightOSM = "d1922b25-af4e-4ba3-84af-fe9bea896051"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MapTiles = "fea52e5a-b371-463b-85f5-81770daa2737"
 OSMMakie = "76b6901f-8821-46bb-9129-841bc9cfe677"
 TileProviders = "263fe934-28e1-4ae9-998a-c2629c5fede6"

docs/src/.vitepress/config.mts

@@ -52,6 +52,7 @@ export default defineConfig({
                         text: "Map3D",
                         link: "/map-3d",
                     },
+                    { text: 'Base maps', link: '/basemap' }
                 ],
             },
 

docs/src/.vitepress/theme/style.css

@@ -234,4 +234,18 @@ kbd {
   border-radius: 6px;
   box-shadow: inset 0 -1px 0 var(--vp-c-bg-soft-mute);
   line-height: 0.95em;
+}
+
+/* Component: Docstring Custom Block */
+
+.jldocstring.custom-block {
+  border: 1px solid var(--vp-c-gray-2);
+  color: var(--vp-c-text-1)
+}
+
+.jldocstring.custom-block summary {
+  font-weight: 700;
+  cursor: pointer;
+  user-select: none;
+  margin: 0 0 8px;
 }

docs/src/basemap.md

@@ -0,0 +1,79 @@
+```@meta
+CollapsedDocStrings = true
+```
+
+# Base maps
+
+A "base map" is essentially a static image composed of map tiles, accessible through the [`basemap`](@ref) function.  
+
+This function returns a tuple of `(x, y, image)`, where `x` and `y` are the dimensions of the image, and `image` is the image itself.  
+The image, and axes, are always in the Web Mercator coordinate system.
+
+While Tyler does not currently work with CairoMakie, this method is a way to add a basemap to a CairoMakie plot, though it will not update on zoom.
+
+```@docs; canonical=false
+basemap
+```
+
+## Examples
+
+
+Below are some cool examples of using basemaps.
+
+#### Cover example of London
+
+````@example coverlondon
+using Tyler, TileProviders
+
+xs, ys, img = basemap(
+    TileProviders.OpenStreetMap(),
+    Extent(X=(-0.5, 0.5), Y=(51.25, 51.75)),
+    (1024, 1024)
+)
+````
+
+````@example coverlondon
+using Makie
+image(xs, ys, img; axis = (; aspect = DataAspect()))
+````
+
+Note that the image is in the Web Mercator projection, as are the axes we see here.
+
+#### NASA GIBS tileset, plotted as a `meshimage` on a `GeoAxis`
+
+````@example nasagibs
+const BACKEND = Makie.current_backend() # hide
+using Tyler, TileProviders, GeoMakie, CairoMakie, Makie
+
+provider = TileProviders.NASAGIBS(:ViirsEarthAtNight2012)
+
+xs, ys, img = basemap(provider, Extent(X=(-90, 90), Y=(-90, 90)), (1024, 1024))
+
+meshimage(
+    xs, ys, img; 
+    source = "+proj=webmerc", # REMEMBER: `img` is always in Web Mercator...
+    axis = (; type = GeoAxis, dest = "+proj=ortho +lat_0=0 +lon_0=0"),
+    npoints = 1024,
+)
+````
+
+````@example nasagibs
+BACKEND.activate!() # hide
+````
+
+
+### OpenSnowMap on polar stereographic projection
+
+````@example opensnowmap
+using Tyler, TileProviders, GeoMakie, Makie
+
+meshimage(
+    basemap(
+        TileProviders.OpenSnowMap(),
+        Extent(X=(-180, 180), Y=(50, 90)),
+        (1024, 1024)
+    )...;
+    source = "+proj=webmerc",
+    axis = (; type = GeoAxis, dest = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-45"),
+)
+````

src/Tyler.jl

@@ -47,6 +47,6 @@ include("tile-fetching.jl")
 include("provider/interpolations.jl")
 include("provider/elevation/elevation-provider.jl")
 include("provider/pointclouds/geotiles-pointcloud-provider.jl")
-
+include("basemap.jl")
 
 end

src/basemap.jl

@@ -0,0 +1,145 @@
+#=
+# Static basemaps
+
+This file provides the ability to get static base maps from Tyler.
+
+Its main entry point is the `basemap` function, which returns a tuple 
+`(x, y, z)` of the image data (`z`) and its axes (`x` and `y`).
+
+This file also contains definitions for `convert_arguments` that make
+the following syntax "just work":
+```julia
+image(TileProviders.Google(), Rect2f(-0.0921, 51.5, 0.04, 0.025), (1000, 1000); axis= (; aspect = DataAspect()))
+```
+
+You do still have to provide the extent and image size, but this is substantially better than nothing.
+
+=#
+
+export basemap
+
+
+"""
+    _z_index(extent::Extent, res::NamedTuple, crs::WebMercator) => Int
+
+Calculate a z value from `extent` and pixel resolution `res` for `crs`.
+The rounded mean calculated z-index for X and Y resolutions is returned.
+
+`res` should be `(X=xres, Y=yres)` to match the extent.
+
+We assume tiles are the standard 256*256 pixels. Note that this is not an
+enforced standard, and that retina tiles are 512*512.
+"""
+function _z_index(extent::Union{Rect,Extent}, res::NamedTuple, crs; tile_size = 256)
+    # Calculate the number of tiles at each z and get the one
+    # closest to the resolution `res`
+    target_ntiles = prod(map(r -> r / tile_size, 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
+
+"""
+    basemap(provider::TileProviders.Provider, bbox::Extent; size, res, min_zoom_level, max_zoom_level)::(xs, ys, img)
+
+Download the most suitable basemap for the given bounding box and size, return a tuple of `(x_interval, y_interval, image)`.
+All returned coordinates and images are in the **Web Mercator** coordinate system (since that's how tiles are defined).
+
+The input bounding box must be in the **WGS84** (long/lat) coordinate system.
+
+## Example
+
+```julia
+basemap(TileProviders.Google(), Extent(X = (-0.0921, -0.0521), Y = (51.5, 51.525)), size   = (1000, 1000))
+```
+
+## Keyword arguments
+
+`size` and `res` are mutually exclusive, and you must provide one.
+
+`min_zoom_level - 0` and `max_zoom_level = 16` are the minimum and maximum zoom levels to consider.
+
+`res` should be a tuple of the form `(X=xres, Y=yres)` to match the extent.
+"""
+function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Rect2{<: Real}, Extent}; 
+        size = nothing, res = nothing, min_zoom_level = 0, max_zoom_level = 16
+    )
+    bbox = Extents.extent(boundingbox)
+    # First, handle keyword arguments
+    @assert (isnothing(size) || isnothing(res)) "You must provide either `size` or `res`, but not both."
+    @assert !(isnothing(size) && isnothing(res)) "You must provide either the `size` or `res` keywords.  Current values: $(size), $(res)"
+    _size = if isnothing(size) 
+        # convert resolution to size using bbox and round(Int, x)
+        (round(Int, (bbox.X[2] - bbox.X[1]) / first(res)), round(Int, (bbox.Y[2] - bbox.Y[1]) / last(res)))
+    else
+        (first(size), last(size))
+    end
+    return basemap(provider, bbox, _size; min_zoom_level, max_zoom_level)
+end
+
+function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Rect2{<: Real}, Extent}, size::Tuple{Int, Int}; min_zoom_level = 0, max_zoom_level = 16)
+    bbox = Extents.extent(boundingbox)
+    # 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.WebMercator())
+    #= TODO:
+    Here we assume all tiles are 256x256.  
+    It's easy to compute this though, by either:
+    - Making a sample query for the tile (0, 0, 0) (but you are not guaranteed this exists)
+    - Some function that returns the tile size for a given provider / dispatch form
+    =#
+    tile_widths = (256, 256)
+    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 = FileIO.load(FileIO.query(IOBuffer(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 also produce the image's axes:
+    # Note that this is in the Web Mercator coordinate system.
+    return (tilegrid_extent.X, tilegrid_extent.Y, image_receptacle)
+end
+
+# We also use this in some Makie converts to allow `image` to work
+Makie.used_attributes(trait::Makie.ImageLike, provider::TileProviders.AbstractProvider, bbox::Union{Rect2, Extent}, size::Union{Int, Tuple{Int, Int}}) = (:min_zoom_level, :max_zoom_level)
+
+function Makie.convert_arguments(trait::Makie.ImageLike, provider::TileProviders.AbstractProvider, bbox::Extent, size::Union{Int, Tuple{Int, Int}}; min_zoom_level = 0, max_zoom_level = 16)
+    return Makie.convert_arguments(trait, basemap(provider, bbox, (first(size), last(size)); min_zoom_level, max_zoom_level)...)
+end
+
+function Makie.convert_arguments(trait::Makie.ImageLike, provider::TileProviders.AbstractProvider, bbox::Rect2, size::Union{Int, Tuple{Int, Int}}; min_zoom_level = 0, max_zoom_level = 16)
+    return Makie.convert_arguments(trait, provider, Extents.extent(bbox), (first(size), last(size)); min_zoom_level, max_zoom_level)
+end
+
+

test/runtests.jl

@@ -38,6 +38,14 @@ display(m)
     @test GeoInterface.crs(m) == Tyler.MapTiles.WebMercator()
 end
 
+@testset "Basemap" begin
+    @test_nowarn Tyler.basemap(Tyler.TileProviders.Google(), london, (1000, 1000))
+    @test_nowarn Tyler.basemap(Tyler.TileProviders.Google(), london; size = (1000, 1000))
+    @test_nowarn Tyler.basemap(Tyler.TileProviders.Google(), london; res = 0.001)
+    x, y, img = Tyler.basemap(Tyler.TileProviders.Google(), london, (1000, 1000))
+    @test img isa Matrix{<: Makie.RGBA}
+end
+
 # Reference tests?
 # provider = TileProviders.NASAGIBS()
 # m = Tyler.Map(Rect2f(0, 50, 40, 20), 5; provider=provider, min_tiles=8, max_tiles=32)