From d22d733b7d436bab904baa1caf28b1e6e7004680 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 21:03:14 +0200 Subject: [PATCH 1/9] Add a function `basemap` and some converts for image --- src/Tyler.jl | 1 + src/basemap.jl | 95 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 8 ++++ 3 files changed, 104 insertions(+) create mode 100644 src/basemap.jl diff --git a/src/Tyler.jl b/src/Tyler.jl index 60dd73a..6d44617 100644 --- a/src/Tyler.jl +++ b/src/Tyler.jl @@ -13,6 +13,7 @@ using OrderedCollections: OrderedCollections, OrderedSet using ThreadSafeDicts: ThreadSafeDicts, ThreadSafeDict using TileProviders: TileProviders, AbstractProvider, geturl, min_zoom, max_zoom +include("basemap.jl") include("interpolations.jl") const TileImage = Matrix{RGB{N0f8}} diff --git a/src/basemap.jl b/src/basemap.jl new file mode 100644 index 0000000..18268db --- /dev/null +++ b/src/basemap.jl @@ -0,0 +1,95 @@ +#= +# 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. + +=# + +""" + basemap(provider::TileProviders.Provider, bbox::Extent; size, res, min_zoom_level = 0, max_zoom_level = 16)::(xs, ys, img) +""" +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." + _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.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 also produce the image's axes: + xs = (..)(tilegrid_extent.X...) + ys = (..)(tilegrid_extent.Y...) + # image(ras; axis = (; aspect = DataAspect())) + return (xs, ys, 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 167be4b..40d0dc7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,6 +38,14 @@ m = wait(Tyler.Map(london; scale=1)) # waits until all tiles are displayed @test GeoInterface.crs(m) == Tyler.MapTiles.WebMercator() end +@testset "Basemap" begin + @test_nowarn basemap(TileProviders.Google(), london, (1000, 1000)) + @test_nowarn basemap(TileProviders.Google(), london; size = (1000, 1000)) + @test_nowarn basemap(TileProviders.Google(), london; res = 0.001) + x, y, img = basemap(TileProviders.Google(), london, (1000, 1000)) + @test img isa Matrix{<: RGBA} +end + # Reference tests? # provider = TileProviders.NASAGIBS() # m = Tyler.Map(Rect2f(0, 50, 40, 20), 5; provider=provider, min_tiles=8, max_tiles=32) From 622cb762443bffe1b473c8547c49435758a501b7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 21:14:39 +0200 Subject: [PATCH 2/9] Make this work properly in a fresh env --- src/basemap.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/basemap.jl b/src/basemap.jl index 18268db..1d3a24b 100644 --- a/src/basemap.jl +++ b/src/basemap.jl @@ -42,6 +42,12 @@ function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Re # Compute the dimensions of the tile grid, so we can feed them into a # Raster later. tilegrid_extent = Extents.extent(tilegrid, MapTiles.WGS84()) + #= 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) + =# + 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. @@ -75,10 +81,7 @@ function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Re end # Now, we have a complete image. # We can also produce the image's axes: - xs = (..)(tilegrid_extent.X...) - ys = (..)(tilegrid_extent.Y...) - # image(ras; axis = (; aspect = DataAspect())) - return (xs, ys, image_receptacle) + return (tilegrid_extent.X, tilegrid_extent.Y, image_receptacle) end # We also use this in some Makie converts to allow `image` to work From 79fd4c5e76aabab3455f84c47fea05ca1999c816 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 21:35:06 +0200 Subject: [PATCH 3/9] Fix the kwarg asserts --- src/basemap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/basemap.jl b/src/basemap.jl index 1d3a24b..59f55e1 100644 --- a/src/basemap.jl +++ b/src/basemap.jl @@ -23,7 +23,7 @@ function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Re 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." + @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))) From 9c373963e136035ad2799165fe30382904f6badc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 21:35:37 +0200 Subject: [PATCH 4/9] Return an image in WebMercator to be consistent Thoughts? How do we determine which CRS the input bbox takes? Its currently WGS84 but do we want to give an option to go to webmerc at all? Conversion between the two is mostly lossless... --- src/basemap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/basemap.jl b/src/basemap.jl index 59f55e1..0886e5f 100644 --- a/src/basemap.jl +++ b/src/basemap.jl @@ -41,7 +41,7 @@ function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Re 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_extent = Extents.extent(tilegrid, MapTiles.WebMercator()) #= TODO: Here we assume all tiles are 256x256. It's easy to compute this though, by either: From 7bde166e298fabb8ce7780386d9ac1618e2f004c Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 22:14:07 +0200 Subject: [PATCH 5/9] Fix tests --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 40d0dc7..423cf94 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,10 +39,10 @@ m = wait(Tyler.Map(london; scale=1)) # waits until all tiles are displayed end @testset "Basemap" begin - @test_nowarn basemap(TileProviders.Google(), london, (1000, 1000)) - @test_nowarn basemap(TileProviders.Google(), london; size = (1000, 1000)) - @test_nowarn basemap(TileProviders.Google(), london; res = 0.001) - x, y, img = basemap(TileProviders.Google(), london, (1000, 1000)) + @test_nowarn basemap(Tyler.TileProviders.Google(), london, (1000, 1000)) + @test_nowarn basemap(Tyler.TileProviders.Google(), london; size = (1000, 1000)) + @test_nowarn basemap(Tyler.TileProviders.Google(), london; res = 0.001) + x, y, img = basemap(Tyler.TileProviders.Google(), london, (1000, 1000)) @test img isa Matrix{<: RGBA} end From e756e9df7e6974a45a4e2293355f3f4bdf738247 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 22:29:47 +0200 Subject: [PATCH 6/9] Really fix tests --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 423cf94..0684f52 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,10 +39,10 @@ m = wait(Tyler.Map(london; scale=1)) # waits until all tiles are displayed end @testset "Basemap" begin - @test_nowarn basemap(Tyler.TileProviders.Google(), london, (1000, 1000)) - @test_nowarn basemap(Tyler.TileProviders.Google(), london; size = (1000, 1000)) - @test_nowarn basemap(Tyler.TileProviders.Google(), london; res = 0.001) - x, y, img = basemap(Tyler.TileProviders.Google(), london, (1000, 1000)) + @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{<: RGBA} end From 2fdfd455024810c6b091749d30652ac9770ad831 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 6 Jul 2024 23:22:07 +0200 Subject: [PATCH 7/9] grrrrrr --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 0684f52..59c195a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -43,7 +43,7 @@ end @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{<: RGBA} + @test img isa Matrix{<: Makie.RGBA} end # Reference tests? From 860040de116c3e1fde356090b672bebe924fca63 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 9 Jul 2024 10:34:18 +0200 Subject: [PATCH 8/9] Use the new FileIO approach --- src/basemap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/basemap.jl b/src/basemap.jl index 0886e5f..7955a6f 100644 --- a/src/basemap.jl +++ b/src/basemap.jl @@ -62,7 +62,7 @@ function basemap(provider::TileProviders.AbstractProvider, boundingbox::Union{Re 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) + img = FileIO.load(FileIO.query(IOBuffer(result))) # 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 = ( From 1a4bcadf894b90909745f8d1352a006b22b71d4b Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 9 Jul 2024 13:23:57 +0200 Subject: [PATCH 9/9] Add a basic `basemap` page (needs a lot of help) --- docs/src/.vitepress/config.mts | 1 + docs/src/basemap.md | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 docs/src/basemap.md diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts index dd41854..1254c54 100644 --- a/docs/src/.vitepress/config.mts +++ b/docs/src/.vitepress/config.mts @@ -45,6 +45,7 @@ export default defineConfig({ { text: 'Whale shark trajectory', link: '/whale_shark' }, { text: 'Ice loss animation', link: '/iceloss_ex' }, { text: 'Interpolation On The Fly', link: '/interpolation' }, + { text: 'Base maps', link: '/basemap' }, ] }, diff --git a/docs/src/basemap.md b/docs/src/basemap.md new file mode 100644 index 0000000..2f76469 --- /dev/null +++ b/docs/src/basemap.md @@ -0,0 +1,10 @@ +# Base maps + +A "base map" is essentially a static image composed of map tiles, accessible through the [`basemap`](@ref) function. + +## Examples + +NASA GIBS tileset on GeoAxis + +OpenSnowMap sounds really cool, maybe a North Pole projection? +