diff --git a/docs/Project.toml b/docs/Project.toml index fd80009..fb62e9c 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts index 8c7a7ff..d7c5b10 100644 --- a/docs/src/.vitepress/config.mts +++ b/docs/src/.vitepress/config.mts @@ -52,6 +52,7 @@ export default defineConfig({ text: "Map3D", link: "/map-3d", }, + { text: 'Base maps', link: '/basemap' } ], }, diff --git a/docs/src/.vitepress/theme/style.css b/docs/src/.vitepress/theme/style.css index 969af9b..7283617 100644 --- a/docs/src/.vitepress/theme/style.css +++ b/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; } \ No newline at end of file diff --git a/docs/src/basemap.md b/docs/src/basemap.md new file mode 100644 index 0000000..b110944 --- /dev/null +++ b/docs/src/basemap.md @@ -0,0 +1,82 @@ +```@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 +using Tyler.Extents + +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 +using Tyler, TileProviders, GeoMakie, CairoMakie, Makie +using Tyler.Extents +const BACKEND = Makie.current_backend() # hide + +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, GLMakie +using Tyler.Extents + +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"), +) +```` diff --git a/src/Tyler.jl b/src/Tyler.jl index 3b3770e..b54c239 100644 --- a/src/Tyler.jl +++ b/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 diff --git a/src/basemap.jl b/src/basemap.jl new file mode 100644 index 0000000..75165f5 --- /dev/null +++ b/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 + + diff --git a/test/runtests.jl b/test/runtests.jl index 58549b5..5754384 100644 --- a/test/runtests.jl +++ b/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)