diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index ed7ab57..d3edf9f 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -21,16 +21,15 @@ permissions: id-token: write statuses: write -# Allow only one concurrent deployment, skipping runs queued between the run in-progress and latest queued. -# However, do NOT cancel in-progress runs as we want to allow these production deployments to complete. +# Cancel previous runs if a new one is triggered (by push on PR) concurrency: - group: pages - cancel-in-progress: false + group: ${{ github.workflow }}-${{ github.head_ref || github.run_id }} + cancel-in-progress: true jobs: # Build job build: - runs-on: ubuntu-latest + runs-on: ubuntu-20.04 steps: - name: Checkout uses: actions/checkout@v4 @@ -51,4 +50,4 @@ jobs: JULIA_DEBUG: "Documenter" DATADEPS_ALWAYS_ACCEPT: true run: | - DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ --color=yes docs/make.jl \ No newline at end of file + DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ --color=yes docs/make.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ba5ab73..bb4ec87 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,7 @@ concurrency: jobs: test: name: Tests - runs-on: ubuntu-latest + runs-on: ubuntu-20.04 env: DISPLAY: ':0' steps: @@ -25,9 +25,7 @@ jobs: version: 1 arch: x64 - uses: julia-actions/cache@v2 - - run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev libcairo2-dev libfreetype6-dev libffi-dev libjpeg-dev libpng-dev libz-dev - - name: Install pkgs dependencies - run: xvfb-run -s '-screen 0 1024x768x24' julia --project=@. -e 'using Pkg; pkg"add https://github.com/JuliaGeo/TileProviders.jl https://github.com/JuliaGeo/MapTiles.jl"' + - run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev xsettingsd x11-xserver-utils - uses: julia-actions/julia-runtest@v1 with: - prefix: xvfb-run -s '-screen 0 1024x768x24' + prefix: xvfb-run -s '-screen 0 1024x768x24' diff --git a/Project.toml b/Project.toml index ecd7d58..dcbc6bf 100644 --- a/Project.toml +++ b/Project.toml @@ -4,16 +4,27 @@ authors = ["SimonDanisch ", "Lazaro Alonso DocumenterVitepress.jl
', + copyright: `© Copyright ${new Date().getUTCFullYear()}. Released under the MIT License.`, + }, }, - socialLinks: [ - { icon: 'github', link: 'https://github.com/MakieOrg/Tyler.jl' } - ], - footer: { - message: 'Made with DocumenterVitepress.jl
', - copyright: `© Copyright ${new Date().getUTCFullYear()}. Released under the MIT License.` - } - } -}) \ No newline at end of file +}); diff --git a/docs/src/assets/alpine.png b/docs/src/assets/alpine.png new file mode 100644 index 0000000..eab2fe5 Binary files /dev/null and b/docs/src/assets/alpine.png differ diff --git a/docs/src/assets/data/logo.png b/docs/src/assets/data/logo.png deleted file mode 100644 index 529a628..0000000 Binary files a/docs/src/assets/data/logo.png and /dev/null differ diff --git a/docs/src/assets/data/favicon.png b/docs/src/assets/favicon.png similarity index 100% rename from docs/src/assets/data/favicon.png rename to docs/src/assets/favicon.png diff --git a/docs/src/assets/data/logo.ico b/docs/src/assets/logo.ico similarity index 100% rename from docs/src/assets/data/logo.ico rename to docs/src/assets/logo.ico diff --git a/docs/src/assets/pointclouds.png b/docs/src/assets/pointclouds.png new file mode 100644 index 0000000..d435249 Binary files /dev/null and b/docs/src/assets/pointclouds.png differ diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index e1753bd..7c33633 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -13,7 +13,7 @@ In the Julia REPL type: ``` The `]` character starts the Julia [package manager](https://docs.julialang.org/en/v1/stdlib/Pkg/). Hit backspace key to return to Julia prompt. -Or, explicitly use `Pkg` +Or, explicitly use `Pkg` ```julia using Pkg @@ -25,12 +25,8 @@ Pkg.add(["Tyler.jl"]) ````@example london using Tyler, GLMakie m = Tyler.Map(Rect2f(-0.0921, 51.5, 0.04, 0.025)) -wait(m) -save("london.png", current_figure()) # hide -nothing # hide ```` -![](london.png) ::: info @@ -46,21 +42,17 @@ We can use a different tile provider as well as any style `theme` from Makie as using GLMakie, Tyler using Tyler.TileProviders -provider = TileProviders.OpenTopoMap() +provider = TileProviders.OpenStreetMap(:Mapnik) london = Rect2f(-0.0921, 51.5, 0.04, 0.025) with_theme(theme_dark()) do m = Tyler.Map(london; provider) hidedecorations!(m.axis) hidespines!(m.axis) - wait(m) + m end -save("londonProvider.png", current_figure()) # hide -nothing # hide ```` -![](londonProvider.png) - ## Providers list More providers are available. See the following list: @@ -91,15 +83,11 @@ with_theme(theme_dark()) do fig = Figure(; size =(1200,600)) ax = Axis(fig[1,1]) # aspect = DataAspect() m = Tyler.Map(london; provider, figure=fig, axis=ax) + wait(m) hidedecorations!(ax) hidespines!(ax) - wait(m) fig end -save("londonFigure.png", current_figure()) # hide -nothing # hide ```` -![](londonFigure.png) - -Next, we could add any other plot type on top of the `ax` axis defined above. \ No newline at end of file +Next, we could add any other plot type on top of the `ax` axis defined above. diff --git a/docs/src/iceloss_ex.md b/docs/src/iceloss_ex.md index 23b4dc9..6190054 100644 --- a/docs/src/iceloss_ex.md +++ b/docs/src/iceloss_ex.md @@ -69,22 +69,15 @@ show map ````@example ice fig = Figure(; size = (1200,600)) ax = Axis(fig[1,1]) -m = Tyler.Map(extent; provider, figure=fig, axis=ax); -wait(m) -save("ice_loss1.png", current_figure()) # hide -nothing # hide +m = Tyler.Map(extent; provider, figure=fig, axis=ax) ```` -![](ice_loss1.png) - create initial scatter plot ````@example ice scatter!(ax, X, Y; color = Z, colormap = cmap, colorrange = [0, n], markersize = 10); -save("ice_loss2.png", current_figure()) # hide -nothing # hide +m ```` -![](ice_loss2.png) add colorbar @@ -98,29 +91,25 @@ Colorbar(fig[1,2]; colormap = cmap, colorrange = [a,b], hidedecorations!(ax); # hide frames hidespines!(ax); -# wait for tiles to fully load -wait(m) -save("ice_loss3.png", current_figure()) # hide -nothing # hide +m ```` -![](ice_loss3.png) -loop to create animation +loop to create animation ````julia -for k = 1:15 +for k = 1:15 # reset apha alpha[:] = zeros(nc); cmap[] = Colors.alphacolor.(cmap[], alpha) - for i in 2:1:n + for i in 2:1:n # modify alpha alpha[1:maximum([1,round(Int64,i*nc/n)])] = alpha[1:maximum([1,round(Int64,i*nc/n)])] .* (1.05^-1.5); alpha[maximum([1,round(Int64,i*nc/n)])] = 1; cmap[] = Colors.alphacolor.(cmap[], alpha); sleep(0.001); - end + end end ```` ```@raw html -``` \ No newline at end of file +``` diff --git a/docs/src/interpolation.md b/docs/src/interpolation.md index 4fd509c..dbb146a 100644 --- a/docs/src/interpolation.md +++ b/docs/src/interpolation.md @@ -32,13 +32,8 @@ p2 = Tyler.Interpolator(fun; options) b = Rect2f(-20.0, -20.0, 40.0, 40.0) m = Tyler.Map(b, provider=p1) -wait(m) -save("interpolation.png", current_figure()) # hide -nothing # hide ```` -![](interpolation.png) - ::: tip Try `b = Rect2f(-180.0, -89.9, 360.0, 179.8)` diff --git a/docs/src/map-3d.md b/docs/src/map-3d.md new file mode 100644 index 0000000..c3271c2 --- /dev/null +++ b/docs/src/map-3d.md @@ -0,0 +1,133 @@ +# Map3D + +Tyler also offers to view tiles in 3D, and offers a simple Elevation and PointCloud provider. + +````@example map3d +using Tyler, GLMakie +using Tyler: ElevationProvider + +lat, lon = (47.087441, 13.377214) +delta = 0.3 +ext = Rect2f(lon-delta/2, lat-delta/2, delta, delta) +m = Tyler.Map3D(ext; provider=ElevationProvider()) +```` + +## Elevation with PlotConfig + +With PlotConfig one can change the way tiles are plotted. +The it has a preprocess + postprocess function and allows to pass any plot attribute to the tile. +These attributes are global and will be passed to every tile plot. + +````@example map3d +lat, lon = (47.087441, 13.377214) +delta = 0.3 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +cfg = Tyler.PlotConfig( + preprocess=pc -> map(p -> p .* 2, pc), + shading=FastShading, + colormap=:alpine +) +m = Tyler.Map3D(ext; provider=ElevationProvider(nothing), plot_config=cfg) +```` + +## PointClouds + +The PointCloud provider downloads from [geotiles.citg.tudelft](https://geotiles.citg.tudelft.nl), which spans most of the netherlands. + +````@example map3d +lat, lon = (52.40459835, 4.84763329) +delta = 0.03 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +provider = Tyler.GeoTilePointCloudProvider() +m = Tyler.Map3D(ext; provider=provider) +```` + + +## Pointclouds with PlotConfig + Meshscatter + +There is also a MeshScatter plot config, which can be used to switch the point cloud plotting from scatter to meshscatter. +This looks better, at a significant slow down. + +````@example map3d +lat, lon = (52.40459835, 4.84763329) +delta = 0.008 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +provider = Tyler.GeoTilePointCloudProvider() +image = ElevationProvider(nothing) +cfg = Tyler.MeshScatterPlotconfig() +m1 = Tyler.Map3D(ext; provider=provider, plot_config=cfg) +cfg = Tyler.PlotConfig(preprocess=pc -> map(p -> p .* 2, pc), shading=FastShading, colormap=:alpine) +m2 = Tyler.Map3D(m1; provider=image, plot_config=cfg) +m1 +```` + + +## Using RPRMakie + +```julia +using RPRMakie, FileIO +function render_rpr(m, name, radiance=1000000) + wait(m) + ax = m.axis + cam = ax.scene.camera_controls + lightpos = Vec3f(cam.lookat[][1], cam.lookat[][2], cam.eyeposition[][3]) + lights = [ + EnvironmentLight(1.5, load(RPR.assetpath("studio026.exr"))), + PointLight(lightpos, RGBf(radiance, radiance * 0.9, radiance * 0.9)) + ] + empty!(ax.scene.lights) + append!(ax.scene.lights, lights) + save("$(name).png", ax.scene; plugin=RPR.Northstar, backend=RPRMakie, iterations=2000) +end +function plastic_material() + return (type=:Uber, reflection_color=Vec4f(1), + reflection_weight=Vec4f(1), reflection_roughness=Vec4f(0.1), + reflection_anisotropy=Vec4f(0), reflection_anisotropy_rotation=Vec4f(0), + reflection_metalness=Vec4f(0), reflection_ior=Vec4f(1.4), refraction_weight=Vec4f(0), + coating_weight=Vec4f(0), sheen_weight=Vec4f(0), emission_weight=Vec3f(0), + transparency=Vec4f(0), reflection_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), + emission_mode=UInt(RPR.RPR_UBER_MATERIAL_EMISSION_MODE_SINGLESIDED), + coating_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), sss_multiscatter=true, + refraction_thin_surface=true) +end +``` + +### Elevation with RPRMakie + +```julia +lat, lon = (47.087441, 13.377214) +delta = 0.5 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +cfg = Tyler.PlotConfig( + preprocess=pc -> map(p -> p .* 2, pc), + shading=FastShading, + material=plastic_material(), + colormap=:alpine +) +m = Tyler.Map3D(ext; + provider=ElevationProvider(nothing), + plot_config=cfg, + max_plots=5 +) +render_rpr(m, "alpine", 10000000) +``` + +![](./assets/alpine.png) + +### PointClouds with RPRMakie + +```julia +lat, lon = (52.40459835229174, 4.84763329882317) +delta = 0.005 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +provider = Tyler.GeoTilePointCloudProvider() +mat = (type=:Microfacet, roughness=0.2, ior=1.390) +cfg = Tyler.MeshScatterPlotconfig(material=mat, markersize=4) +m = Tyler.Map3D(ext; provider=provider, plot_config=cfg, max_plots=3, size=(2000, 2000)) +cfg = Tyler.PlotConfig(material=mat, colormap=:Blues) +m2 = Tyler.Map3D(m; provider=ElevationProvider(nothing), plot_config=cfg, max_plots=5) +render_rpr(m, "pointclouds") +cp(Tyler) +``` + +![](./assets/pointclouds.png) diff --git a/docs/src/osmmakie.md b/docs/src/osmmakie.md index aa9a085..f07e9ae 100644 --- a/docs/src/osmmakie.md +++ b/docs/src/osmmakie.md @@ -1,4 +1,4 @@ -## OpenStreetMap data (OSM) +## OpenStreetMap data (OSM) In this example, we combine OpenStreetMap data, loading some roads and buildings and plotting them on top of a Tyler map. @@ -38,6 +38,5 @@ m = Tyler.Map(london; provider=provider, crs=Tyler.wgs84) m.axis.aspect = map_aspect(area.minlat, area.maxlat) p = osmplot!(m.axis, osm; buildings) # DataInspector(m.axis) # this is broken/slow -wait(m) - +m ```` diff --git a/docs/src/points_poly_text.md b/docs/src/points_poly_text.md index 404613e..323d754 100644 --- a/docs/src/points_poly_text.md +++ b/docs/src/points_poly_text.md @@ -32,31 +32,25 @@ set how much area to map in degrees and define an `Extent` for display in web_me ````@example plottypes delta = 1 -extent = Extent(X = (lon - delta/2, lon + delta/2), Y = (lat-delta/2, lat+delta/2)); +extent = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta); ```` show map ````@example plottypes -m = Tyler.Map(extent; provider, figure=Figure(; size=(1000, 600))) -# wait for tiles to fully load -wait(m) -save("map_plottypes.png", current_figure()) # hide -nothing # hide +m = Tyler.Map(extent; provider, size=(1000, 600)) ```` -![](map_plottypes.png) - now plot a point, polygon and text on the map ````@example plottypes objscatter = scatter!(m.axis, pts; color = :red, marker = '⭐', markersize = 50) # hide ticks, grid and lables -hidedecorations!(m.axis) +hidedecorations!(m.axis) # hide frames hidespines!(m.axis) -# Plot a plygon on the map +# Plot a plygon on the map p1 = (lon-delta/8, lat-delta/8) p2 = (lon-delta/8, lat+delta/8) p3 = (lon+delta/8, lat+delta/8) @@ -72,7 +66,4 @@ text!(pts2, text = "Basic Example"; fontsize = 30, color = :darkblue, align = (:center, :center) ) m -save("map_plottypes_all.png", current_figure()) # hide -nothing # hide ```` -![](map_plottypes_all.png) diff --git a/docs/src/whale_shark.md b/docs/src/whale_shark.md index 862cc9d..8973119 100644 --- a/docs/src/whale_shark.md +++ b/docs/src/whale_shark.md @@ -21,6 +21,7 @@ function to_web_mercator(lo,lat) end nothing # hide ```` + and downloading and preparing the data is done next ````@example whale @@ -56,11 +57,7 @@ fig = Figure(; size = (1200, 600)) ax = Axis(fig[1,1]) m = Tyler.Map(Rect2f(Rect2f(lomn - δlon/2, lamn-δlat/2, 2δlon, 2δlat)); provider, figure=fig, axis=ax) -wait(m) -save("whale_init.png", fig) # hide -nothing # hide ```` -![](whale_init.png) ## Initial point @@ -79,11 +76,9 @@ objscatter = scatter!(ax, whale; markersize = 15, color = :orangered, strokecolor=:white, strokewidth=1.5) hidedecorations!(ax) hidespines!(ax) - -save("whale_init_point.png", fig) # hide -nothing # hide +m ```` -![](whale_init_point.png) + ## Animated trajectory @@ -104,4 +99,4 @@ set_theme!() # reset theme (default) ```@raw html -``` \ No newline at end of file +``` diff --git a/src/3d-map.jl b/src/3d-map.jl new file mode 100644 index 0000000..18d1509 --- /dev/null +++ b/src/3d-map.jl @@ -0,0 +1,141 @@ + +using GeometryBasics, LinearAlgebra + +function setup_axis!(axis::LScene, ext_target, crs) + # Disconnect all events + scene = axis.scene + cam = scene.camera_controls + Makie.disconnect!(scene.camera) + # add back our mouse controls + add_tyler_mouse_controls!(scene, cam) + X = ext_target.X + Y = ext_target.Y + pmin = Vec3(X[1], Y[1], 0) + xrange = X[2] - X[1] + center = Vec3((X[1] + X[2]) / 2, (Y[1] + Y[2]) / 2, 0) + xaxis = Vec3(xrange / 2, 0, 0) + eyeposition = pmin .+ xaxis .+ Vec3(0, 0, 0.8 * xrange) + up = cross(xaxis, center .- eyeposition) + cam.settings.clipping_mode[] = :static + update_cam!(axis.scene, eyeposition, center, up) + return +end + +function Map3D(extent, extent_crs=wgs84; + size=(1000, 1000), + figure=Makie.Figure(; size=size), + axis=Makie.LScene(figure[1, 1]; show_axis=false), + provider=TileProviders.OpenStreetMap(:Mapnik), + fetching_scheme=Tiling3D(), + args... + ) + return Map( + extent, extent_crs; + figure=figure, + axis=axis, + provider=provider, + fetching_scheme=fetching_scheme, + args... + ) +end + +function tile_reloader(m::Map{LScene}) + scene = m.axis.scene + camc = scene.camera_controls + update_signal = map((a,b)->nothing, camc.lookat, camc.eyeposition) + throttled = Makie.Observables.throttle(0.2, update_signal) + on(scene, throttled; update=true) do _ + update_tiles!(m, (scene.camera, camc)) + return + end +end + + +function frustum_snapshot(cam::Camera) + bottom_left = Point4d(-1, -1, 1, 1) + top_left = Point4d(-1, 1, 1, 1) + top_right = Point4d(1, 1, 1, 1) + bottom_right = Point4d(1, -1, 1, 1) + rect_ps = [bottom_left, top_left, top_right, bottom_right] + inv_pv = inv(cam.projectionview[]) + return map(rect_ps) do p + p = inv_pv * p + return p[Vec(1, 2, 3)] / p[4] + end +end + +# Function to find the intersection of a ray with a plane +function ray_plane_intersection(position::Point3, normal::Vec3, start::Point3, direction::Vec3)::Union{Point3,Nothing} + # Calculate the denominator of the intersection formula + denom = dot(normal, direction) + if abs(denom) > 1e-6 + # Calculate the numerator of the intersection formula + numerator = dot(normal, position .- start) + t = numerator / denom + if t >= 0 + # The intersection point is at t * direction from the ray's start point + intersection = start .+ t .* direction + return intersection + end + end + # Return nothing if there is no intersection + return nothing +end + +function area_around_lookat(camc::Camera3D, multiplier=2) + x, y, _ = camc.lookat[] + z = camc.eyeposition[][3] + w = multiplier * z + return Rect2d(x - w, y - w, 2w, 2w) +end + +function frustrum_plane_intersection(cam::Camera, camc::Camera3D) + frustrum = frustum_snapshot(cam) + result = Union{Nothing, Point3d}[] + eyepos = camc.eyeposition[] + for point in frustrum + res = ray_plane_intersection(Point3d(0), Vec3d(0, 0, 1), Point3d(eyepos), Vec3d(point .- eyepos)) + push!(result, res) + end + x, y, _ = camc.lookat[] + z = camc.eyeposition[][3] + w = 2 * z + fallback = [Point3d(x - w, y - w, 0), Point3d(x - w, y + w, 0), Point3d(x + w, y + w, 0), Point3d(x + w, y - w, 0)] + return map(enumerate(result)) do (i, res) + if !isnothing(res) + return res + else + return fallback[i] + end + end +end + +function to_rect(extent::Extent) + return Rect2(extent.X[1], extent.Y[1], extent.X[2] - extent.X[1], extent.Y[2] - extent.Y[1]) +end + +function tiles_from_poly(m::Map, points) + mini = minimum(points) + points_s = sort(points, by=(x -> norm(x .- mini))) + start = points_s[1] + diagonal = norm(points_s[end] .- start) + zoom, actual_zoom, ntiles_approx = optimal_zoom(m, diagonal) + if abs(actual_zoom - zoom) > 2 || ntiles_approx > m.max_plots + @warn "Too many tiles to plot, which means zoom level is not supported. Plotting no tiles for this zoomlevel." maxlog = 1 + return OrderedSet{Tile}() + end + m.zoom[] = actual_zoom + boundingbox = Rect2d(Rect3d(points)) + ext = Extents.extent(boundingbox) + tilegrid = TileGrid(ext, actual_zoom, m.crs) + tiles = OrderedSet{Tile}() + poly = Polygon(map(p-> Point2d(p[1], p[2]), points)) + for tile in tilegrid + tile_ext = to_rect(Extents.extent(tile, m.crs)) + tile_poly = Polygon(decompose(Point2d, tile_ext)) + if GO.intersects(poly, tile_poly) + push!(tiles, tile) + end + end + return tiles +end diff --git a/src/Tyler.jl b/src/Tyler.jl index 60dd73a..382390d 100644 --- a/src/Tyler.jl +++ b/src/Tyler.jl @@ -1,439 +1,52 @@ module Tyler -using Colors: Colors, RGB, N0f8 +using Colors: Colors, RGB, N0f8, Colorant using Extents: Extents, Extent using GeoInterface: GeoInterface -using GeometryBasics: GeometryBasics, GLTriangleFace, Point2f, Vec2f, Rect2f, Rect2, Rect, decompose, decompose_uv +using GeometryBasics: GeometryBasics, GLTriangleFace, Point2f, Vec2f, Rect2f, Rect2, Rect, decompose, decompose_uv, Vec3d, Point2d, Point3d, Point4d using HTTP: HTTP using LRUCache: LRUCache, LRU using MapTiles: MapTiles, Tile, TileGrid, web_mercator, wgs84, CoordinateReferenceSystemFormat -using Makie: Makie, Observable, Figure, Axis, RGBAf, on, isopen, linesegments!, meta, mesh!, translate!, scale! -using Makie: FileIO, ImageIO +using Makie: Makie, Observable, Figure, Axis, LScene, RGBAf, on, isopen, meta, mesh!, translate!, scale!, Plot using OrderedCollections: OrderedCollections, OrderedSet using ThreadSafeDicts: ThreadSafeDicts, ThreadSafeDict using TileProviders: TileProviders, AbstractProvider, geturl, min_zoom, max_zoom +using Makie +using Makie: AbstractAxis +using LinearAlgebra, GeometryBasics +using GeometryBasics +using Proj +using Statistics, DelimitedFiles +using PointClouds +using ArchGDAL +import GeoFormatTypes as GFT +using Downloads +using Scratch +using ImageIO, FileIO +import GeometryOps as GO + + +const CACHE_PATH = Ref("") + +function __init__() + # Initialize at init for relocatability + CACHE_PATH[] = @get_scratch!("download-cache") +end + +abstract type AbstractPlotConfig end +abstract type FetchingScheme end +abstract type AbstractMap end + +include("downloader.jl") +include("tiles.jl") +include("map.jl") +include("3d-map.jl") +include("tyler-cam3d.jl") +include("tile-plotting.jl") +include("tile-fetching.jl") +include("provider/interpolations.jl") +include("provider/elevation/elevation-provider.jl") +include("provider/pointclouds/geotiles-pointcloud-provider.jl") -include("interpolations.jl") - -const TileImage = Matrix{RGB{N0f8}} - -""" - Map - - Map(extent, [extent_crs=wgs84]; kw...) - -Tylers main object, it plots tiles onto a Makie.jl `Axis`, -downloading and plotting more tiles as you zoom and pan. - -# Arguments - --`extent`: the initial extent of the map, as a `GeometryBasics.Rect` - or an `Extents.Extent` in the projection of `extent_crs`. --`extent_crs`: Any `GeoFormatTypes` compatible crs, the default is wsg84. - -# Keywords - --`resolution`: The figure resolution. --`figure`: an existing `Makie.Figure` object. --`crs`: The providers coordinate reference system. --`provider`: a TileProviders.jl `Provider`. --`max_parallel_downloads`: limits the attempted simultaneous downloads, with a default of `16`. --`cache_size_gb`: limits the cache for storing tiles, with a default of `5`. --`depth`: the number of layers to load when zooming. Lower numbers will be slightly faster - but have more artefacts. The default is `8`. --`halo`: The fraction of the width of tiles to add as a halo so that panning is smooth - the - tiles will already be loaded. The default is `0.2`, which means `0.1` on each side. --`scale`: a tile scaling factor. Low number decrease the downloads but reduce the resolution. - The default is `1.0`. -""" -struct Map{A<:Makie.AbstractAxis} - provider::AbstractProvider - crs::CoordinateReferenceSystemFormat - zoom::Observable{Int} - figure::Figure - axis::A - displayed_tiles::OrderedSet{Tile} - plots::Dict{Tile,Any} - free_tiles::Vector{Makie.Plot} - fetched_tiles::LRU{Tile,TileImage} - max_parallel_downloads::Int - # TODO, use Channel here - queued_but_not_downloaded::OrderedSet{Tile} - tiles_being_added::ThreadSafeDict{Tile,Task} - downloaded_tiles::Channel{Tuple{Tile,TileImage}} - display_task::Base.RefValue{Task} - download_task::Base.RefValue{Task} - depth::Int - halo::Float64 - scale::Float64 - max_zoom::Int -end - -# Wait for all tiles to be loaded -function Base.wait(map::Map) - # The download + plot loops need a screen to do their work! - if isnothing(Makie.getscreen(map.figure.scene)) - display(map.figure) - end - screen = Makie.getscreen(map.figure.scene) - isnothing(screen) && error("No screen after display. Wrong backend? Only WGLMakie and GLMakie are supported.") - while true - if !isempty(map.tiles_being_added) - wait(last(first(map.tiles_being_added))) - end - if !isempty(map.queued_but_not_downloaded) - sleep(0.001) # we don't have a task to wait on, so we sleep - end - # We're done if both are empty! - if isempty(map.tiles_being_added) && isempty(map.queued_but_not_downloaded) - while !isempty(map.downloaded_tiles) - sleep(0.01) - end - return map - end - end -end - -Base.showable(::MIME"image/png", ::Map) = true -function Base.show(io::IO, m::MIME"image/png", map::Map) - wait(map) - Makie.show(io, m, map.figure) -end - -function stopped_displaying(screen::Makie.MakieScreen) - Backend = parentmodule(typeof(screen)) - if nameof(Backend) == :WGLMakie - session = Backend.get_screen_session(screen) - isnothing(session) && return false - !isready(session) && return false - return !isopen(session) - elseif nameof(Backend) == :GLMakie - return !isopen(screen) - else - error("Unsupported backend: $(Backend)") - end -end - -function stopped_displaying(fig::Figure) - scene = Makie.get_scene(fig) - screen = Makie.getscreen(scene) - # if not displayed yet, we return true, since we're using - # is_open as a condition in our while loop to stop once it got displayed & closed - isnothing(screen) && return false - return stopped_displaying(screen) -end - -function Map(extent, extent_crs=wgs84; - resolution=(1000, 1000), - figure=Makie.Figure(; size=resolution), - axis=Makie.Axis(figure[1, 1]; aspect=Makie.DataAspect()), - provider=TileProviders.OpenStreetMap(:Mapnik), - crs=MapTiles.web_mercator, - max_parallel_downloads=16, - cache_size_gb=5, - depth=8, - halo=0.2, - scale=2.0, - max_zoom=TileProviders.max_zoom(provider) -) - - fetched_tiles = LRU{Tile, Matrix{RGB{N0f8}}}(; maxsize=cache_size_gb * 10^9, by=Base.sizeof) - free_tiles = Makie.Combined[] - tiles_being_added = ThreadSafeDict{Tile,Task}() - downloaded_tiles = Channel{Tuple{Tile,TileImage}}(128) - - display_task = Base.RefValue{Task}() - download_task = Base.RefValue{Task}() - - # Extent - # if extent input is a HyperRectangle then convert to type Extent - extent isa Extent || (extent = Extents.extent(extent)) - ext_target = MapTiles.project_extent(extent, extent_crs, crs) - X = ext_target.X - Y = ext_target.Y - if axis isa Makie.Axis3 - Makie.xlims!(axis, X[1], X[2]) - Makie.ylims!(axis, Y[1], Y[2]) - elseif axis isa Makie.LScene - # do nothing - else # any other axis is 2d - axis.autolimitaspect = 1 - Makie.limits!(axis, (X[1], X[2]), (Y[1], Y[2])) - end - plots = Dict{Tile,Any}() - tyler = Map( - provider, crs, - Observable(1), - figure, axis, OrderedSet{Tile}(), plots, free_tiles, - fetched_tiles, - max_parallel_downloads, OrderedSet{Tile}(), - tiles_being_added, downloaded_tiles, - display_task, download_task, - depth, Float64(halo), Float64(scale), max_zoom - ) - tyler.zoom[] = get_zoom(tyler, extent) - download_task[] = @async begin - while !stopped_displaying(figure) - # we dont download all tiles at once, so when one download task finishes, we may want to schedule more downloads: - if !isempty(tyler.queued_but_not_downloaded) - queue_tile!(tyler, popfirst!(tyler.queued_but_not_downloaded)) - end - sleep(0.01) - end - empty!(tyler.queued_but_not_downloaded) - empty!(tyler.displayed_tiles) - @debug("stopped download task") - end - display_task[] = @async begin - while !stopped_displaying(figure) - while isnothing(Makie.getscreen(tyler.axis.scene)) - sleep(0.01) - end - screen = Makie.getscreen(tyler.axis.scene) - while !isopen(screen) - sleep(0.01) - end - tile, img = take!(downloaded_tiles) - try - create_tile_plot!(tyler, tile, img) - # fixes on demand renderloop which doesn't pick up all updates! - if hasproperty(screen, :requires_update) - screen.requires_update = true - end - catch e - @warn "error while creating tile" exception = (e, Base.catch_backtrace()) - end - end - @debug("stopped display task") - end - - # Queue tiles to be downloaded & displayed - update_tiles!(tyler, ext_target) - - if axis isa Union{Makie.Axis, Makie.Axis3} - on(axis.scene, axis.finallimits) do extent - stopped_displaying(figure) && return - update_tiles!(tyler, extent) - return - end - else - # we can't hook it up so don't do anything. - # the user can zoom manually using `update_tiles!`. - end - return tyler -end - -GeoInterface.crs(tyler::Map) = tyler.crs -Extents.extent(tyler::Map) = Extents.extent(Makie.Rect2d(tyler.axis.finallimits[])) - -function stop_download!(map::Map, tile::Tile) - # delete!(map.tiles_being_added, tile) - # TODO can we actually interrupt downloads? - # Doesn't seem to work this way at least: - # if haskey(map.tiles_being_added, tile) - # task = map.tiles_being_added[tile] - # ex = InterruptException() - # Base.throwto(task, ex) - # end -end - -function remove_tiles!(map::Map, tiles_being_displayed::OrderedSet{Tile}) - to_remove_plots = setdiff(keys(map.plots), tiles_being_displayed) - for tile in to_remove_plots - if haskey(map.plots, tile) - plot = pop!(map.plots, tile) - plot.visible = false - push!(map.free_tiles, plot) - end - end - remove_from_queue = setdiff(map.queued_but_not_downloaded, tiles_being_displayed) - foreach(t-> delete!(map.queued_but_not_downloaded, t), remove_from_queue) - to_remove_downloads = setdiff(keys(map.tiles_being_added), tiles_being_displayed) - foreach(t -> stop_download!(map, t), to_remove_downloads) -end - -function create_tileplot!(axis, image) - # Use a mesh plot until image rotr90 is fixed - # Also, this will make it easier to subdivide the image mesh - # To apply some other coordinate transforms - # use GeometryBasics.Tesselation(rect, (128, 128)) to get a 128x128 subdivied mesh - rect = Rect2f(0, 0, 1, 1) - points = decompose(Point2f, rect) - faces = decompose(GLTriangleFace, rect) - uv = decompose_uv(rect) - map!(uv -> Vec2f(uv[1], 1 - uv[2]), uv, uv) - m = GeometryBasics.Mesh(meta(points; uv=uv), faces) - # Plot directly into scene to not update limits - return mesh!(axis.scene, m; color=image, shading=Makie.NoShading, inspectable=false) -end - -function place_tile!(tile::Tile, plot, crs) - bounds = MapTiles.extent(tile, crs) - xmin, xmax = bounds.X - ymin, ymax = bounds.Y - translate!(plot, xmin, ymin, tile.z - 100) - scale!(plot, xmax - xmin, ymax - ymin, 1) - return -end - -function create_tile_plot!(tyler::Map, tile::Tile, image::TileImage) - if haskey(tyler.plots, tile) - # this shouldn't get called with plots that are already displayed - @debug "getting tile plot already plotted" - return tyler.plots[tile] - end - if isempty(tyler.free_tiles) - mplot = create_tileplot!(tyler.axis, image) - else - mplot = pop!(tyler.free_tiles) - mplot.color[] = image - mplot.visible = true - end - tyler.plots[tile] = mplot - place_tile!(tile, mplot, tyler.crs) -end - -function fetch_tile(tyler::Map, tile::Tile) - #println("Fetching Tile") - return get!(tyler.fetched_tiles, tile) do - fetch_tile(tyler.provider, tile) - end -end -function fetch_tile(provider::AbstractProvider, tile::Tile) - url = TileProviders.geturl(provider, tile.x, tile.y, tile.z) - result = HTTP.get(url; retry=false, readtimeout=4, connect_timeout=4) - io = IOBuffer(result.body) # this wraps the byte data in an IO-like type - format = FileIO.query(io) # this interrogates the magic bits to see what file format it is (JPEG, PNG, etc) - return FileIO.load(format) # this actually loads the data using ImageIO.jl or whatever other FileIO loader exists -end - -function queue_tile!(tyler::Map, tile) - queue = tyler.tiles_being_added - # NO need to start a need task! - haskey(queue, tile) && return - if length(queue) > tyler.max_parallel_downloads - # queue them for being downloaded once all other downloads are finished - # This helps not spamming the server and also gives a chance to delete - # the download request when e.g. zooming out really fast - push!(tyler.queued_but_not_downloaded, tile) - else - queue[tile] = @async try - # TODO, i think we should better check if we want to display this tile, - # even though we just scheduled it - if (tile in tyler.displayed_tiles) - img = fetch_tile(tyler, tile) - # we may have moved already and the tile doesn't need to be displayed anymore - if tile in tyler.displayed_tiles - put!(tyler.downloaded_tiles, (tile, img)) - end - end - catch e - if !(e isa InterruptException) - @warn "error while downloading tile" exception=e - end - # if the tile still needs to be displayed, reque it - # TODO should only retry for certain errors - if tile in tyler.displayed_tiles - @info("requeing after download error!") - delete!(queue, tile) - queue_tile!(tyler, tile) - end - finally - delete!(queue, tile) - end - end -end - -TileProviders.max_zoom(tyler::Map) = tyler.max_zoom -TileProviders.min_zoom(tyler::Map) = Int(min_zoom(tyler.provider)) - -function get_zoom(tyler::Map, area) - res = size(tyler.axis.scene) .* tyler.scale - clamp(z_index(area, (X=res[2], Y=res[1]), tyler.crs), min_zoom(tyler), max_zoom(tyler)) -end - -function update_tiles!(tyler::Map, area::Union{Rect,Extent}) - area = typeof(area) <: Rect ? Extents.extent(area) : area - # `depth` determines the number of layers below the current - # layer to load. Tiles are downloaded in order from lowest to highest zoom. - depth = tyler.depth - - # Calculate the zoom level - zoom = get_zoom(tyler, area) - tyler.zoom[] = zoom - - # And the z layers we will plot - layer_range = max(min_zoom(tyler), zoom - depth):zoom - # Get the tiles around the mouse first - xpos, ypos = Makie.mouseposition(tyler.axis.scene) - xspan = (area.X[2] - area.X[1]) * 0.01 - yspan = (area.Y[2] - area.Y[1]) * 0.01 - mouse_area = Extents.Extent(X=(xpos - xspan, xpos + xspan), Y=(ypos - yspan, ypos + yspan)) - # Make a halo around the mouse tile to load next, intersecting area so we don't download outside the plot - mouse_halo_area = grow_extent(mouse_area, 10) - # Define a halo around the area to download last, so pan/zoom are filled already - halo_area = grow_extent(area, tyler.halo) # We don't mind that the middle tiles are the same, the OrderedSet will remove them - # Define all the tiles in the order they will load in - areas = if Extents.intersects(mouse_halo_area, area) - mha = Extents.intersection(mouse_halo_area, area) - if Extents.intersects(mouse_area, area) - [Extents.intersection(mouse_area, area), mha, area, halo_area] - else - [mha, area, halo_area] - end - else - [area, halo_area] - end - - area_layers = map(layer_range) do z - map(areas) do ext - MapTiles.TileGrid(ext, z, tyler.crs) - end |> Iterators.flatten - end |> Iterators.flatten - - # Create the full set of tiles - # Using an ordered set gives a smoother load than a Set - new_tiles_set = OrderedSet{Tile}(area_layers) - # Remove any tiles not in the new set - remove_tiles!(tyler, new_tiles_set) - - to_add = setdiff(new_tiles_set, tyler.displayed_tiles) - - # replace - empty!(tyler.displayed_tiles) - union!(tyler.displayed_tiles, new_tiles_set) - - # Queue tiles to be downloaded & displayed - foreach(tile -> queue_tile!(tyler, tile), to_add) -end - -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 grow_extent(area::Union{Rect,Extent}, factor) - map(Extents.bounds(area)) do axis_bounds - span = axis_bounds[2] - axis_bounds[1] - pad = factor * span / 2 - (axis_bounds[1] - pad, axis_bounds[2] + pad) - end |> Extent -end - -function debug_tile!(map::Tyler.Map, tile::Tile) - plot = Makie.linesegments!(map.axis, Rect2f(0, 0, 1, 1), color=:red, linewidth=1) - Tyler.place_tile!(tile, plot, web_mercator) -end - -function debug_tiles!(map::Tyler.Map) - for tile in map.displayed_tiles - debug_tile!(map, tile) - end -end end diff --git a/src/downloader.jl b/src/downloader.jl new file mode 100644 index 0000000..bcea7c5 --- /dev/null +++ b/src/downloader.jl @@ -0,0 +1,54 @@ + +abstract type AbstractDownloader end + +struct NoDownload <: AbstractDownloader end + +struct ByteDownloader <: AbstractDownloader + timeout::Float64 + downloader::Downloads.Downloader + io::IOBuffer + bytes::Vector{UInt8} +end + +function ByteDownloader(timeout=3) + downloader = Downloads.Downloader() + downloader.easy_hook = (easy, info) -> Downloads.Curl.setopt(easy, Downloads.Curl.CURLOPT_LOW_SPEED_TIME, timeout) + return ByteDownloader(timeout, downloader, IOBuffer(), UInt8[]) +end + +function download_tile_data(dl::ByteDownloader, provider, url) + Downloads.download(url, dl.io; downloader=dl.downloader, timeout=dl.timeout) + # a bit of shananigans to allocate less and stress the GC less! + resize!(dl.bytes, dl.io.ptr - 1) + copyto!(dl.bytes, 1, dl.io.data, 1, dl.io.ptr-1) + seekstart(dl.io) + return dl.bytes +end + +struct PathDownloader <: AbstractDownloader + timeout::Float64 + downloader::Downloads.Downloader + cache_dir::String + lru::LRU{String, Int} +end + +function PathDownloader(cache_dir; timeout=5, cache_size_gb=5) + isdir(cache_dir) || mkpath(cache_dir) + lru = LRU{String, Int}(maxsize=cache_size_gb * 10^9, by=identity) + downloader = Downloads.Downloader() + downloader.easy_hook = (easy, info) -> Downloads.Curl.setopt(easy, Downloads.Curl.CURLOPT_LOW_SPEED_TIME, timeout) + return PathDownloader(timeout, downloader, cache_dir, lru) +end + +function unique_filename(url) + return string(hash(url)) +end + +function download_tile_data(dl::PathDownloader, provider::AbstractProvider, url) + unique_name = unique_filename(url) + path = joinpath(dl.cache_dir, unique_name * file_ending(provider)) + if !isfile(path) + Downloads.download(url, path; downloader=dl.downloader) + end + return path +end diff --git a/src/map.jl b/src/map.jl new file mode 100644 index 0000000..0eb2b34 --- /dev/null +++ b/src/map.jl @@ -0,0 +1,267 @@ + +function tile_key(provider::AbstractProvider, tile::Tile) + return TileProviders.geturl(provider, tile.x, tile.y, tile.z) +end + + +""" + Map(extent, [extent_crs=wgs84]; kw...) + Map(map::Map; ...) # layering another provider on top of an existing map + +Tylers main object, it plots tiles onto a Makie.jl `Axis`, +downloading and plotting more tiles as you zoom and pan. +When layering providers over each other with `Map(map::Map; ...)`, you can use `toggle_visibility!(map)` to hide/unhide them. + +# Arguments + +-`extent`: the initial extent of the map, as a `GeometryBasics.Rect` + or an `Extents.Extent` in the projection of `extent_crs`. +-`extent_crs`: Any `GeoFormatTypes` compatible crs, the default is wsg84. + +# Keywords + +-`size`: The figure size. +-`figure`: an existing `Makie.Figure` object. +-`crs`: The providers coordinate reference system. +-`provider`: a TileProviders.jl `Provider`. +-`max_parallel_downloads`: limits the attempted simultaneous downloads, with a default of `16`. +-`cache_size_gb`: limits the cache for storing tiles, with a default of `5`. +-`fetching_scheme=Halo2DTiling()`: The tile fetching scheme. Can be SimpleTiling(), Halo2DTiling(), or Tiling3D(). +-`scale`: a tile scaling factor. Low number decrease the downloads but reduce the resolution. + The default is `0.5`. +-`plot_config`: A `PlotConfig` object to change the way tiles are plotted. +-`max_zoom`: The maximum zoom level to display, with a default of `TileProviders.max_zoom(provider)`. +-`max_plots=400:` The maximum number of plots to keep displayed at the same time. +""" +struct Map{Ax<:Makie.AbstractAxis} <: AbstractMap + provider::AbstractProvider + figure::Figure + axis::Ax + plot_config::AbstractPlotConfig + # The tile downloader + cacher + tiles::TileCache + # The tiles for the current zoom level - we may plot many more than this + current_tiles::ThreadSafeDict{Tile,Bool} + # The plots we have created but are not currently visible and can be reused + unused_plots::Vector{Makie.Plot} + # All tile plots we're currently plotting + plots::ThreadSafeDict{String,Tuple{Makie.Plot,Tile,Rect}} + # All tiles we currently wish to be plotting, but may not yet be downloaded + displayed + should_get_plotted::ThreadSafeDict{String,Tile} + display_task::Base.RefValue{Task} + + crs::CoordinateReferenceSystemFormat + zoom::Observable{Int} + fetching_scheme::FetchingScheme + max_zoom::Int + max_plots::Int + scale::Float64 +end + +setup_axis!(::Makie.AbstractAxis, ext_target, crs) = nothing + +function setup_axis!(axis::Axis, ext_target, crs) + X = ext_target.X + Y = ext_target.Y + axis.autolimitaspect = 1 + Makie.limits!(axis, (X[1], X[2]), (Y[1], Y[2])) + axis.elements[:background].depth_shift[] = 0.1f0 + translate!(axis.elements[:background], 0, 0, -1000) + axis.elements[:background].color = :transparent + axis.xgridvisible = false + axis.ygridvisible = false + return +end + +function Map3D(m::Map; kw...) + ax = m.axis + # Make a copy of the lscene, so we can easier separate the plots. + ax2 = LScene(ax.parent, ax.layoutobservables, ax.blockscene) + ax2.scene = Scene(ax.scene; camera=ax.scene.camera, camera_controls=ax.scene.camera_controls) + return Map3D(nothing, nothing; figure=m.figure, axis=ax2, kw...) +end + + +""" + Map(m::Map; kw...) +Layering constructor to show another provider on top of an existing map. + +## Example +```julia +lat, lon = (52.395593, 4.884704) +delta = 0.01 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +m1 = Tyler.Map(ext) +m2 = Tyler.Map(m1; provider=TileProviders.Esri(:WorldImagery), plot_config=Tyler.PlotConfig(alpha=0.5, postprocess=(p-> translate!(p, 0, 0, 1f0)))) +m1 +``` +""" +function Map(m::Map; kw...) + ax = m.axis + # Make a copy of the lscene, so we can easier separate the plots. + ax2 = Axis(ax.parent, ax.layoutobservables, ax.blockscene) + ax2.scene = Scene(ax.scene; camera=ax.scene.camera, camera_controls=ax.scene.camera_controls) + setfield!(ax2, :elements, ax.elements) + setfield!(ax2, :targetlimits, ax.targetlimits) + setfield!(ax2, :finallimits, ax.finallimits) + setfield!(ax2, :block_limit_linking, ax.block_limit_linking) + setfield!(ax2.scene, :float32convert, ax.scene.float32convert) + return Map(nothing, nothing; figure=m.figure, axis=ax2, kw...) +end + +toggle_visibility!(m::Map) = m.axis.scene.visible[] = !m.axis.scene.visible[] + +function Base.close(m::Map) + cleanup_queue!(m, OrderedSet{Tile}()) + empty!(m.current_tiles) + empty!(m.unused_plots) + empty!(m.plots) + empty!(m.should_get_plotted) + close(m.tiles) +end + +function Map(extent, extent_crs=wgs84; + size=(1000, 1000), + figure=Makie.Figure(; size=size), + axis=Makie.Axis(figure[1, 1]; aspect=Makie.DataAspect()), + plot_config=PlotConfig(), + provider=TileProviders.OpenStreetMap(:Mapnik), + crs=MapTiles.web_mercator, + cache_size_gb=5, + max_parallel_downloads=1, + fetching_scheme=Halo2DTiling(), + max_zoom=TileProviders.max_zoom(provider), + max_plots=400, + scale=1) + + # Extent + # if extent input is a HyperRectangle then convert to type Extent + ext_target = nothing + if !isnothing(extent) && !isnothing(extent_crs) + extent isa Extent || (extent = Extents.extent(extent)) + ext_target = MapTiles.project_extent(extent, extent_crs, crs) + setup_axis!(axis, ext_target, crs) + end + + tiles = TileCache(provider; cache_size_gb=cache_size_gb, max_parallel_downloads=max_parallel_downloads) + downloaded_tiles = tiles.downloaded_tiles + + plots = ThreadSafeDict{String,Tuple{Makie.Plot,Tile,Rect}}() + should_get_plotted = ThreadSafeDict{String,Tile}() + current_tiles = ThreadSafeDict{Tile,Bool}() + unused_plots = Makie.Plot[] + display_task = Base.RefValue{Task}() + + map = Map( + provider, figure, + axis, + plot_config, + tiles, + current_tiles, + unused_plots, + plots, + should_get_plotted, + display_task, crs, + Observable(1), + fetching_scheme, max_zoom, max_plots, Float64(scale) + ) + + map.zoom[] = 0 + closed = Threads.Atomic{Bool}(false) + + display_task[] = @async for (tile, data) in downloaded_tiles + closed[] && break + try + if isnothing(data) + # download went wrong or provider doesn't have tile. + # That means we won't plot this tile and it should not be in the queue anymore + delete!(map.should_get_plotted, tile_key(map.provider, tile)) + else + create_tile_plot!(map, tile, data) + end + catch e + @warn "error while creating tile" exception = (e, Base.catch_backtrace()) + end + end + + tile_reloader(map) + + on(axis.scene.events.window_open) do open + if !open && !closed[] + closed[] = true + # remove all queued tiles! + cleanup_queue!(map, OrderedSet{Tile}()) + close(map) + end + end + return map +end + +function Base.wait(m::AbstractMap; timeout=50) + # The download + plot loops need a screen to do their work! + if isnothing(Makie.getscreen(m.figure.scene)) + screen = display(m.figure.scene) + end + screen = Makie.getscreen(m.figure.scene) + isnothing(screen) && error("No screen after display.") + wait(m.tiles; timeout=timeout) + start = time() + while true + tile_keys = Set(tile_key.((m.provider,), keys(m.current_tiles))) + if all(k -> haskey(m.plots, k), tile_keys) + break + end + if time() - start > timeout + @warn "Timeout waiting for all tiles to be plotted" + break + end + sleep(0.01) + end + return m +end + +function tile_reloader(map::Map{Axis}) + axis = map.axis + throttled = Makie.Observables.throttle(0.2, axis.finallimits) + on(axis.scene, throttled; update=true) do extent + update_tiles!(map, extent) + return + end +end + +function plotted_tiles(m::Map) + return getindex.(values(m.plots), 2) +end + +GeoInterface.crs(map::Map) = map.crs +Extents.extent(map::Map) = Extents.extent(map.axis.finallimits[]) + +TileProviders.max_zoom(map::Map) = map.max_zoom +TileProviders.min_zoom(map::Map) = Int(min_zoom(map.provider)) + +Base.showable(::MIME"image/png", ::AbstractMap) = true + +function Base.show(io::IO, m::MIME"image/png", map::AbstractMap) + wait(map) + Makie.show(io, m, map.figure.scene) +end + +function Base.display(map::AbstractMap) + wait(map) + Base.display(map.figure.scene) +end + + +function remove_unused!(m::AbstractMap, tile::Tile) + return remove_unused!(m, tile_key(m.provider, tile)) +end + +function remove_unused!(m::AbstractMap, key::String) + plot_tile = get(m.plots, key, nothing) + if !isnothing(plot_tile) + plot, tile, bounds = plot_tile + move_to_back!(plot, abs(m.zoom[] - tile.z), bounds) + return plot, key + end + return nothing +end diff --git a/src/provider/elevation/elevation-provider.jl b/src/provider/elevation/elevation-provider.jl new file mode 100644 index 0000000..e359ed6 --- /dev/null +++ b/src/provider/elevation/elevation-provider.jl @@ -0,0 +1,55 @@ +""" + ElevationProvider(color_provider::Union{Nothing, AbstractProvider}=TileProviders.Esri(:WorldImagery); cache_size_gb=5) + +Provider rendering elevation data from [arcgis](https://elevation3d.arcgis.com/arcgis/rest/services/WorldElevation3D/Terrain3D/ImageServer). +This provider is special, since it uses a second provider for color information, +which also means you can provide a cache size, since color tile caching has to be managed by the provider. +When set to `nothing`, no color provider is used and the elevation data is used to color the surface with a colormap directly. +Use `Map(..., plot_config=Tyler.PlotConfig(colormap=colormap))` to set the colormap and other `surface` plot attributes. +""" +struct ElevationProvider <: AbstractProvider + color_provider::Union{Nothing, AbstractProvider} + tile_cache::LRU{String} + downloader::Vector +end + +Tyler.get_tile_format(::ElevationProvider) = Tyler.ElevationData +TileProviders.options(::ElevationProvider) = nothing +TileProviders.min_zoom(::ElevationProvider) = 0 +TileProviders.max_zoom(::ElevationProvider) = 16 + +function ElevationProvider(provider=TileProviders.Esri(:WorldImagery); cache_size_gb=5) + TileFormat = get_tile_format(provider) + fetched_tiles = LRU{String,TileFormat}(; maxsize=cache_size_gb * 10^9, by=Base.sizeof) + downloader = [get_downloader(provider) for i in 1:Threads.maxthreadid()] + ElevationProvider(provider, fetched_tiles, downloader) +end + +function TileProviders.geturl(::ElevationProvider, x::Integer, y::Integer, z::Integer) + return "https://elevation3d.arcgis.com/arcgis/rest/services/WorldElevation3D/Terrain3D/ImageServer/tile/$(z)/$(y)/$(x)" +end + +function get_downloader(::ElevationProvider) + cache_dir = joinpath(CACHE_PATH[], "ElevationProvider") + return PathDownloader(cache_dir) +end + +file_ending(::ElevationProvider) = ".lerc" + +function fetch_tile(provider::ElevationProvider, dl::PathDownloader, tile::Tile) + path = download_tile_data(dl, provider, TileProviders.geturl(provider, tile.x, tile.y, tile.z)) + dataset = ArchGDAL.read(path, options=["DATATYPE=Float32"]) + band = ArchGDAL.getband(dataset, 1) + mini = -450 + maxi = 8700 + elevation_img = collect(reverse(band; dims=2)) + elevation_img .= Float32.(elevation_img) # .* (maxi - mini) .+ mini + if isnothing(provider.color_provider) + return Tyler.ElevationData(elevation_img, Matrix{RGBf}(undef, 0, 0), Vec2d(mini, maxi)) + end + foto_img = get!(provider.tile_cache, path) do + dl = provider.downloader[Threads.threadid()] + fetch_tile(provider.color_provider, dl, tile) + end + return Tyler.ElevationData(elevation_img, rotr90(foto_img), Vec2d(mini, maxi)) +end diff --git a/src/interpolations.jl b/src/provider/interpolations.jl similarity index 82% rename from src/interpolations.jl rename to src/provider/interpolations.jl index 2881fe8..6d07c47 100644 --- a/src/interpolations.jl +++ b/src/provider/interpolations.jl @@ -14,16 +14,23 @@ struct Interpolator{F} <: AbstractProvider colormap::Vector{RGBAf} options::Dict end + Interpolator(f; colormap=:thermal, options=Dict(:minzoom=>1, :maxzoom=>19)) = Interpolator(f, Makie.to_colormap(colormap), options) -function fetch_tile(interpolator::Interpolator, tile::Tile) +function TileProviders.geturl(::Interpolator, x::Integer, y::Integer, z::Integer) + return "$x,$y,$z" +end + +get_downloader(::Interpolator) = NoDownload() + +function fetch_tile(interpolator::Interpolator, ::NoDownload, tile::Tile) (lon, lat) = _tile2positions(tile) z = permutedims(interpolator.interpolator.(lon, lat)) return [_col(interpolator, i) for i in z] end -# TODO just use Makie plotting for colors, +# TODO just use Makie plotting for colors, # we just need to pass the args throught to it _col(i::Interpolator, x) = RGBAf(Makie.interpolated_getindex(i.colormap, x)) _col(::Interpolator, x::RGBAf) = x @@ -35,7 +42,8 @@ _tile2lng(x, z) = (x / 2^z * 360) - 180 _tile2lat(y, z) = -180 / pi * atan(0.5 * (exp(pi - 2 * pi * y / 2^z) - exp(2 * pi * y / 2^z - pi))) _tile2positions(tile::Tile) = _tile2positions(tile.x, tile.y, tile.z) -function _tile2positions(x, y, z) + +function _tile2positions(x, y, z) rng = range(0.5 / 232, 231.5 / 232,232) lons = [_tile2lng(x + i, z) for i in rng, j in rng] lats = [_tile2lat(y + j, z) for i in rng, j in rng] diff --git a/src/provider/pointclouds/generate-mapping-refined.jl b/src/provider/pointclouds/generate-mapping-refined.jl new file mode 100644 index 0000000..ff8df7c --- /dev/null +++ b/src/provider/pointclouds/generate-mapping-refined.jl @@ -0,0 +1,88 @@ +using DelimitedFiles +using GeometryBasics +using GeoInterface, TileProviders +import GeoDataFrames as GDF +using GeoDataFrames +import GeoFormatTypes as GFT +using Statistics, LinearAlgebra +using Makie: Point2d +using ThreadSafeDicts: ThreadSafeDict +using MapTiles +using Extents, Tyler +import GeometryOps as GO + +p = Point2f[(0, 0), (1, 0), (1, 1), (0, 1)] +GO.intersects(Polygon(p), Polygon(map(x -> x .- 1.1, p))) + +function generate_mapping(web_tile_polys, geom_polys, n_neibhbours, geom_df) + mapping = ThreadSafeDict{Tuple{Int,Int},String}() + Threads.@threads for (i, poly) in tuple.(1:length(web_tile_polys), web_tile_polys) + nclose_idx = n_neibhbours[i] + isempty(nclose_idx) && continue + GO.difference(Polygon(p), Polygon(map(x -> x .- 0.5, p)); target=PolygonTrait()) |> GO.area + npoly_indices = findall(x -> GO.intersects(x, poly), geom_polys[nclose_idx]) + poly_indices = nclose_idx[npoly_indices] + if length(poly_indices) > 1 + sort!(poly_indices; by= x -> GO.area(GO.intersection(poly,geom_polys[x])), rev=true) + end + tile = tilegrid[i] + idx = poly_indices[1] + mapping[(tile.x, tile.y)] = geom_df.GT_AHNSUB[idx] + end + return mapping +end + +function get_ntiles(middle, midpoints, n) + findall(x -> norm(x .- middle) < n, midpoints) +end +path = joinpath(@__DIR__, "netherlands", "AHN_subunits_GeoTiles", "AHN_subunits_GeoTiles.shp") +df = GDF.read(path) +geometry = reproject(df.geometry, GFT.EPSG(28992), GFT.EPSG(3857)) +geo_polys = [Polygon(Point2d.(GeoInterface.coordinates(p)[1])) for p in geometry] +mid_points = [mean(coordinates(p)) for p in geo_polys] +bounds = mapreduce(x -> Rect2f(decompose(Point2f, x)), Base.union, geo_polys) + +tilegrid = TileGrid(Extents.extent(bounds), 16, MapTiles.web_mercator) +web_tile_polys = map(tilegrid) do tile + tile_ext = Tyler.to_rect(MapTiles.extent(tile, Tyler.web_mercator)) + return Polygon(decompose(Point2d, tile_ext)[[1, 3, 4, 2]]) +end +first_poly = geo_polys[1] +rect = Rect2f(coordinates(first_poly)) +n_neighbours = norm(widths(rect)) * 2 +all_n_neighbours = map(web_tile_polys) do poly + nclose_idx = get_ntiles(mean(coordinates(poly)), mid_points, n_neighbours) + return nclose_idx +end + +mapping = generate_mapping(web_tile_polys, geo_polys, all_n_neighbours, df) + +mkeys = collect(keys(mapping)) +mvalues = collect(values(mapping)) +k1 = first.(mkeys) +k2 = last.(mkeys) +matrix = hcat(k1, k2) + +open(joinpath(@__DIR__, "mapping.bin"), "w") do io + write(io, length(mkeys)) + write(io, UInt16.(matrix)) + for s in mkeys + write(io, s) + write(io, UInt8(0)) + end +end + +function read_mapping(path) + open(path, "r") do io + n = read(io, Int) + keys = Matrix{UInt16}(undef, n, 2) + read!(io, keys) + strings = Vector{String}(undef, n) + for i in 1:n + strings[i] = readuntil(io, Char(0)) + end + return Dict(Tuple(Int.(keys[i, :])) => strings[i] for i in 1:n) + end +end + +read_mapping(joinpath(@__DIR__, "mapping.bin")) diff --git a/src/provider/pointclouds/generate-mapping.jl b/src/provider/pointclouds/generate-mapping.jl new file mode 100644 index 0000000..ab35a7b --- /dev/null +++ b/src/provider/pointclouds/generate-mapping.jl @@ -0,0 +1,58 @@ +using DelimitedFiles +using GeometryBasics +using GeoInterface, TileProviders +import GeoDataFrames as GDF +using GeoDataFrames +import GeoFormatTypes as GFT +using Statistics, LinearAlgebra +using Makie: Point2d +using ThreadSafeDicts: ThreadSafeDict +using MapTiles +using Extents +import GeometryOps as OP + +function generate_mapping(tilegrid, mid_points, geom_df) + mapping = ThreadSafeDict{Tuple{Int,Int},String}() + Threads.@threads for tile in tilegrid + tb = MapTiles.extent(tile, Tyler.web_mercator) + wx = (tb.X[2] - tb.X[1]) + wy = (tb.Y[2] - tb.Y[1]) + tile_mid = Point2d(tb.X[1] + wx / 2, tb.Y[1] + wy / 2) + val, idx = findmin(x -> norm(tile_mid .- x), mid_points) + if val < max(wx, wy) + mapping[(tile.x, tile.y)] = geom_df.GT_AHNSUB[idx] + end + end + return mapping +end + +function get_tiles(area, point_rects, tiles) + tilenames = String[] + for (i, points) in enumerate(point_rects) + if any(p -> p in area, points) + push!(tilenames, tiles.GT_AHNSUB[i]) + end + end + return tilenames +end + +path = joinpath(@__DIR__, "netherlands", "AHN_subunits_GeoTiles", "AHN_subunits_GeoTiles.shp") +df = GDF.read(path) +geometry = reproject(df.geometry, GFT.EPSG(28992), GFT.EPSG(3857)) +mid_points = [mean(Point2d.(GeoInterface.coordinates(p)[1])) for p in geometry] +bounds = Rect2d(mid_points) +tilegrid = TileGrid(Extents.extent(bounds), 15, MapTiles.web_mercator) + +mapping = generate_mapping(tilegrid, mid_points, df) + +mkeys = collect(keys(mapping)) +mvalues = collect(values(mapping)) +k1 = first.(mkeys) +k2 = last.(mkeys) +matrix = hcat(k1, k2, mvalues) + +DelimitedFiles.writedlm("mapping.csv", matrix, ',') + +matrix = DelimitedFiles.readdlm("mapping.csv", ',') +mapping2 = Dict((r[1], r[2]) => r[3] for r in eachrow(matrix)) +@assrt mapping2 == mapping diff --git a/src/provider/pointclouds/geotiles-pointcloud-provider.jl b/src/provider/pointclouds/geotiles-pointcloud-provider.jl new file mode 100644 index 0000000..4dc2d78 --- /dev/null +++ b/src/provider/pointclouds/geotiles-pointcloud-provider.jl @@ -0,0 +1,96 @@ + +""" + GeoTilePointCloudProvider(subset="AHN1_T") + +The PointCloud provider downloads from [geotiles.citg.tudelft](https://geotiles.citg.tudelft.nl), which spans most of the netherlands. +You can specify the subset to download from, which can be one of the following: +- AHN1_T (default): The most corse dataset, but also the fastest to download (1-5mb compressed per tile) +- AHN2_T: More detailed dataset (~70mb per tile) +- AHN3_T: ~250mb per tile +- AHN4_T: 300-500mb showing much detail, takes a long time to load each tile (over 1 minute per tile). Use `max_plots=5` to limit the number of tiles loaded at once. +""" +struct GeoTilePointCloudProvider <: TileProviders.AbstractProvider + baseurl::String + subset::String + lookup + projs::Vector{Proj.Transformation} +end + +get_tile_format(::GeoTilePointCloudProvider) = PointCloudData +TileProviders.options(::GeoTilePointCloudProvider) = nothing +TileProviders.min_zoom(::GeoTilePointCloudProvider) = 16 +TileProviders.max_zoom(::GeoTilePointCloudProvider) = 16 + +const AHN_SUB_MAPPING = Dict{Tuple{Int,Int},String}() + +function read_mapping(path) + open(path, "r") do io + n = read(io, Int) + keys = Matrix{UInt16}(undef, n, 2) + read!(io, keys) + strings = Vector{String}(undef, n) + for i in 1:n + strings[i] = readuntil(io, Char(0)) + end + return Dict(Tuple(Int.(keys[i, :])) => strings[i] for i in 1:n) + end +end + +function get_ahn_sub_mapping() + if isempty(AHN_SUB_MAPPING) + path = joinpath(@__DIR__, "mapping.bin") + dict = read_mapping(path) + merge!(AHN_SUB_MAPPING, dict) + end + return AHN_SUB_MAPPING +end + +function GeoTilePointCloudProvider(; baseurl="https://geotiles.citg.tudelft.nl", subset="AHN1_T") + projs = [Proj.Transformation(GFT.EPSG(28992), GFT.EPSG(3857)) for i in 1:Threads.maxthreadid()] + return GeoTilePointCloudProvider(baseurl, subset, get_ahn_sub_mapping(), projs) +end + +function get_downloader(::GeoTilePointCloudProvider) + cache_dir = joinpath(CACHE_PATH[], "GeoTilePointCloudProvider") + return PathDownloader(cache_dir) +end + +file_ending(::GeoTilePointCloudProvider) = ".laz" + +function TileProviders.geturl(p::GeoTilePointCloudProvider, x::Integer, y::Integer, z::Integer) + if z == TileProviders.min_zoom(p) && haskey(p.lookup, (x, y)) + return string(p.baseurl, "/", p.subset, "/", p.lookup[(x, y)], ".LAZ") + end + return nothing +end + +function get_points(las, offset, scale, proj) + return map(las.points) do p + point = Point3f(offset .+ Point3f(p.coords) .* scale) + Point3f(proj(point.data)) + end +end + +function load_tile_data(provider::GeoTilePointCloudProvider, path::String) + pc = LAS(path) + isempty(pc.points) && return nothing + proj = provider.projs[Threads.threadid()] + points = get_points(pc, Point3f(pc.coord_offset), Point3f(pc.coord_scale), proj) + extrema = Point3f.(proj.(Point3f[pc.coord_min, pc.coord_max])) + if hasproperty(pc.points[1], :color_channels) + color = map(pc.points) do p + c = map(x -> N0f8(x / 255), p.color_channels) + return RGB(c[1], c[2], c[3]) + end + else + color = last.(points) # z as fallback + end + best_markersize = Dict( + "AHN1_T" => 9.0, + "AHN2_T" => 5.0, + "AHN3_T" => 4.0, + "AHN4_T" => 2.0 + ) + bb = Rect3d(extrema[1], extrema[2] .- extrema[1]) + return PointCloudData(points, color, bb, best_markersize[provider.subset]) +end diff --git a/src/provider/pointclouds/mapping.bin b/src/provider/pointclouds/mapping.bin new file mode 100644 index 0000000..b76de22 Binary files /dev/null and b/src/provider/pointclouds/mapping.bin differ diff --git a/src/tile-fetching.jl b/src/tile-fetching.jl new file mode 100644 index 0000000..d129bd4 --- /dev/null +++ b/src/tile-fetching.jl @@ -0,0 +1,237 @@ + +function queue_plot!(m::Map, tile) + key = tile_key(m.provider, tile) + # Provider doesn't have a tile for this + isnothing(key) && return + m.should_get_plotted[key] = tile + put!(m.tiles.tile_queue, tile) + return +end + +function cleanup_queue!(m::AbstractMap, to_keep::OrderedSet{Tile}) + queue = m.tiles.tile_queue + lock(queue) do + tiles = Tile[] + queued = queue.data + filter!(queued) do tile + if !(tile in to_keep) + Base._increment_n_avail(queue, -1) + push!(tiles, tile) + return false + else + return true + end + end + for tile in tiles + delete!(m.should_get_plotted, tile_key(m.provider, tile)) + end + end +end + +function update_tiles!(m::Map, arealike) + # Get the tiles to be plotted from the fetching scheme and arealike + new_tiles_set, background_tiles = get_tiles_for_area(m, m.fetching_scheme, arealike) + if length(new_tiles_set) > m.max_plots + @warn "Too many tiles to plot, which means zoom level is not supported. Plotting no tiles for this zoomlevel." maxlog = 1 + new_tiles_set = OrderedSet{Tile}() + background_tiles = OrderedSet{Tile}() + end + queued_or_plotted = values(m.should_get_plotted) + to_add = setdiff(new_tiles_set, queued_or_plotted) + + # We don't add any background tile to the current_tiles, so they stay shifted to the back + # They get added async, so at this point `to_add` won't be in currently_plotted yet + will_be_plotted = union(new_tiles_set, queued_or_plotted) + # replace + empty!(m.current_tiles) + for tile in new_tiles_set + m.current_tiles[tile] = true + end + + # Move all plots to the back, that aren't in the newest tileset anymore + for (key, (plot, tile, bounds)) in m.plots + dist = abs(m.zoom[] - tile.z) + if haskey(m.current_tiles, tile) + move_in_front!(plot, dist, bounds) + else + move_to_back!(plot, dist, bounds) + end + end + + # Queue tiles to be downloaded & displayed + to_add_background = setdiff(background_tiles, will_be_plotted) + # Remove any item from queue, that isn't in the new set + to_keep = union(background_tiles, will_be_plotted) + # Remove all tiles that are not in the new set from the queue + cleanup_queue!(m, to_keep) + + # The unique is needed to avoid tiles referencing the same tile + # TODO, we should really consider to disallow this for tile providers, + # This is currently only allowed because of the PointCloudProvider + background = unique(t -> tile_key(m.provider, t), to_add_background) + foreground = unique(t -> tile_key(m.provider, t), to_add) + + # We lock the queue, to put all tiles in one go into the tile queue + # Since download workers take the last tiles first, foreground tiles go last + # Without the lock, a few (n_download_threads) background tiles would be downloaded first, + # since they will be the last in the queue until we add the foreground tiles + lock(m.tiles.tile_queue) do + foreach(tile -> queue_plot!(m, tile), background) + foreach(tile -> queue_plot!(m, tile), foreground) + end +end + +######################################################################################### +##### Halo2DTiling + +struct Halo2DTiling <: FetchingScheme + depth::Int + halo::Float64 + pixel_scale::Float64 +end + +Halo2DTiling(; depth=8, halo=0.2, pixel_scale=2.0) = Halo2DTiling(depth, halo, pixel_scale) + +function get_tiles_for_area(m::Map{Axis}, scheme::Halo2DTiling, area::Union{Rect,Extent}) + area = typeof(area) <: Rect ? Extents.extent(area) : area + # `depth` determines the number of layers below the current + # layer to load. Tiles are downloaded in order from lowest to highest zoom. + depth = scheme.depth + + # Calculate the zoom level + # TODO, also early return if too many tiles to plot? + ideal_zoom, zoom, approx_ntiles = optimal_zoom(m, norm(widths(to_rect(area)))) + m.zoom[] = zoom + + # And the z layers we will plot + layer_range = max(min_zoom(m), zoom - depth):zoom + # Get the tiles around the mouse first + xpos, ypos = Makie.mouseposition(m.axis.scene) + xspan = (area.X[2] - area.X[1]) * 0.01 + yspan = (area.Y[2] - area.Y[1]) * 0.01 + mouse_area = Extents.Extent(; X=(xpos - xspan, xpos + xspan), Y=(ypos - yspan, ypos + yspan)) + # Make a halo around the mouse tile to load next, intersecting area so we don't download outside the plot + mouse_halo_area = grow_extent(mouse_area, 10) + # Define a halo around the area to download last, so pan/zoom are filled already + halo_area = grow_extent(area, scheme.halo) # We don't mind that the middle tiles are the same, the OrderedSet will remove them + # Define all the tiles in the order they will load in + background_areas = if Extents.intersects(mouse_halo_area, area) + mha = Extents.intersection(mouse_halo_area, area) + if Extents.intersects(mouse_area, area) + [Extents.intersection(mouse_area, area), mha, halo_area] + else + [mha, halo_area] + end + else + [halo_area] + end + + foreground = OrderedSet{Tile}(MapTiles.TileGrid(area, zoom, m.crs)) + background = OrderedSet{Tile}() + for z in layer_range + z == zoom && continue + for ext in background_areas + union!(background, MapTiles.TileGrid(ext, z, m.crs)) + end + end + return foreground, background +end + +######################################################################################### +##### SimpleTiling + +struct SimpleTiling <: FetchingScheme +end + +function get_tiles_for_area(m::Map, ::SimpleTiling, area::Union{Rect,Extent}) + area = typeof(area) <: Rect ? Extents.extent(area) : area + diag = norm(widths(to_rect(area))) + # Calculate the zoom level + ideal_zoom, zoom, approx_ntiles = optimal_zoom(m, diag) + m.zoom[] = zoom + return OrderedSet{Tile}(MapTiles.TileGrid(area, zoom, m.crs)), OrderedSet{Tile}() +end + +######################################################################################### +##### Tiling3D + +struct Tiling3D <: FetchingScheme +end + +function get_tiles_for_area(m::Map{LScene}, ::Tiling3D, (cam, camc)::Tuple{Camera,Camera3D}) + points = frustrum_plane_intersection(cam, camc) + eyepos = camc.eyeposition[] + maxdist, _ = findmax(p -> norm(p[3] .- eyepos), points) + camc.far[] = maxdist + camc.near[] = eyepos[3] * 0.01 + update_cam!(m.axis.scene) + return tiles_from_poly(m, points), OrderedSet{Tile}() +end + +function get_tiles_for_area(m::Map{LScene}, s::SimpleTiling, (cam, camc)::Tuple{Camera,Camera3D}) + area = area_around_lookat(camc) + return get_tiles_for_area(m, s, area) +end + + +######################################################################################### +##### Helper functions + +function get_resolution(map::Map) + screen = Makie.getscreen(map.axis.scene) + return isnothing(screen) ? size(map.axis.scene) .* 1.5 : size(screen) +end + +function grow_extent(area::Union{Rect,Extent}, factor) + Extent(map(Extents.bounds(area)) do axis_bounds + span = axis_bounds[2] - axis_bounds[1] + pad = factor * span / 2 + return (axis_bounds[1] - pad, axis_bounds[2] + pad) + end) +end + +function optimal_zoom(m::Map, diagonal) + diagonal_res = norm(get_resolution(m)) * m.scale + # Go over complete known zoomrange of any provider. + # So that we can get the theoretical optimal zoom level, even if the provider doesn't support it, + # which we can then use to calculate the distance to the supported zoomlevel and may decide to not plot anything. + # (TODO, how exactly can we get this over all providers?) + zoomrange = 1:22 + z = optimal_zoom(m.crs, diagonal, diagonal_res, zoomrange, m.zoom[]) + actual_zoom = clamp(z, min_zoom(m), max_zoom(m)) + return z, actual_zoom, approx_tiles(m, actual_zoom, diagonal) +end + +function optimal_zoom(crs, diagonal, diagonal_resolution, zoom_range, old_zoom) + # TODO, this should come from provider + tile_diag_res = norm((255, 255)) + target_ntiles = diagonal_resolution / tile_diag_res + canditates_dict = Dict{Int,Float64}() + candidates = @NamedTuple{z::Int, ntiles::Float64}[] + for z in zoom_range + ext = Extents.extent(Tile(0, 0, z), crs) + mini, maxi = Point2.(ext.X, ext.Y) + diag = norm(maxi .- mini) + ntiles = diagonal / diag + canditates_dict[z] = ntiles + push!(candidates, (; z, ntiles)) + end + if haskey(canditates_dict, old_zoom) # for the first invokation, old_zoom is 0, which is not a candidate + old_ntiles = canditates_dict[old_zoom] + # If the old zoom level is close to the target number of tiles, return it + # to change the zoom level less often + if old_ntiles > (target_ntiles - 1) && old_ntiles < (target_ntiles + 1) + return old_zoom + end + end + dist, idx = findmin(x -> abs(x.ntiles - target_ntiles), candidates) + return candidates[idx].z +end + +function approx_tiles(m::Map, zoom, diagonal) + ext = Extents.extent(Tile(0, 0, zoom), m.crs) + mini, maxi = Point2.(ext.X, ext.Y) + diag = norm(maxi .- mini) + ntiles_diag = diagonal / diag + return (ntiles_diag / sqrt(2)) ^ 2 +end diff --git a/src/tile-plotting.jl b/src/tile-plotting.jl new file mode 100644 index 0000000..708e1a1 --- /dev/null +++ b/src/tile-plotting.jl @@ -0,0 +1,349 @@ +# don't shift 3d plots +move_in_front!(plot, amount, ::Rect3) = nothing +function move_in_front!(plot, amount, ::Rect2) + if !hasproperty(plot, :depth_shift) + translate!(plot, 0, 0, amount) + else + plot.depth_shift = -amount / 100f0 + end +end + +move_to_back!(plot, amount, ::Rect3) = nothing +function move_to_back!(plot, amount, ::Rect2) + if !hasproperty(plot, :depth_shift) + translate!(plot, 0, 0, -amount) + else + plot.depth_shift = amount / 100f0 + end +end + +struct PlotConfig <: AbstractPlotConfig + attributes::Dict{Symbol, Any} + preprocess::Function + postprocess::Function +end + +""" + PlotConfig(; preprocess=identity, postprocess=identity, plot_attributes...) + +Creates a `PlotConfig` object to influence how tiles are being plotted. +* preprocess(tile_data): Function to preprocess the data before plotting. For a tile provider returning image data, preprocess will be called on the image data before plotting. +* postprocess(tile_data): Function to mutate the plot object after creation. Can be used like this: `(plot)-> translate!(plot, 0, 0, 1)`. +* plot_attributes: Additional attributes to pass to the plot + +## Example + +```julia +using Tyler, GLMakie + +config = PlotConfig( + preprocess = (data) -> data .+ 1, + postprocess = (plot) -> translate!(plot, 0, 0, 1), + color = :red +) +lat, lon = (52.395593, 4.884704) +delta = 0.1 +extent = Extent(; X=(lon - delta / 2, lon + delta / 2), Y=(lat - delta / 2, lat + delta / 2)) +Tyler.Map(extent; provider=Tyler.TileProviders.Esri(:WorldImagery), plot_config=config) +``` +""" +PlotConfig(; preprocess=identity, postprocess=identity, plot_attributes...) = PlotConfig(Dict{Symbol,Any}(plot_attributes), preprocess, postprocess) + + +function get_bounds(tile::Tile, data, crs) + bounds = MapTiles.extent(tile, crs) + return to_rect(bounds) +end + +function remove_plot!(m::Map, key::String) + if !haskey(m.plots, key) + @warn "deleting non-existing plot" + delete!(m.should_get_plotted, key) + return + end + plot, _, _ = m.plots[key] + plot.visible = false + delete!(m.plots, key) + delete!(m.should_get_plotted, key) + push!(m.unused_plots, plot) +end + +get_preprocess(config) = identity +get_preprocess(config::PlotConfig) = config.preprocess +get_postprocess(config) = identity +get_postprocess(config::PlotConfig) = config.postprocess + + +function filter_overlapping!(m::Map, bounds::Rect2, tile, key) + return false + # dont filter for 2d plots +end + +function filter_overlapping!(m::Map, bounds::Rect3, tile, key) + # for 3d meshes, we need to remove any plot in the same 2d area + bounds2d = Rect2d(bounds) + for (other_key, (plot, other_tile, other_bounds)) in copy(m.plots) + other_bounds2d = Rect2d(other_bounds) + # If overlap + if bounds2d in other_bounds2d || other_bounds2d in bounds2d + if haskey(m.current_tiles, tile) + # the new plot has priority since it's in the newest current tile set + remove_plot!(m, other_key) + elseif haskey(m.current_tiles, other_tile) + delete!(m.should_get_plotted, key) + # the existing plot has priority so we skip the new plot + return true + else + # If both are not in current_tiles, we remove the plot farthest away from the current zoom level + if abs(tile.z - m.zoom[]) <= abs(other_tile.z - m.zoom[]) + remove_plot!(m, other_key) + else + delete!(m.should_get_plotted, key) + return true + end + end + end + end + return false +end + +function cull_plots!(m::Map) + if length(m.plots) >= (m.max_plots - 1) + # remove the oldest plot + p_tiles = plotted_tiles(m) + available_to_remove = setdiff(p_tiles, keys(m.current_tiles)) + sort!(available_to_remove, by=tile -> abs(tile.z - m.zoom[])) + n_avail = length(available_to_remove) + need_to_remove = min(n_avail, length(m.plots) - m.max_plots) + to_remove = available_to_remove[1:need_to_remove] + for tile in to_remove + plot_key = remove_unused!(m, tile) + if !isnothing(plot_key) + remove_plot!(m, plot_key[2]) + end + end + end +end + +function create_tile_plot!(m::AbstractMap, tile::Tile, data) + # For providers which have to map the same data to different tiles + # Or providers that have e.g. additional parameters like a date + # the url is a much better key than the tile itself + # TODO, should we instead have custom tiles for certain provider? + key = tile_key(m.provider, tile) + # This can happen for tile providers with overlapping data that doesn't map 1:1 to tiles + if haskey(m.plots, key) + delete!(m.should_get_plotted, key) + return + end + + cfg = m.plot_config + data_processed = get_preprocess(cfg)(data) + bounds = get_bounds(tile, data_processed, m.crs) + + this_got_filtered = filter_overlapping!(m, bounds, tile, key) + this_got_filtered && return # skip plotting if it overlaps with a more important plot + + # Cull plots over plot limit + cull_plots!(m) + + if isempty(m.unused_plots) + mplot = create_tileplot!(cfg, m.axis, data_processed, bounds, (tile, m.crs)) + else + mplot = pop!(m.unused_plots) + update_tile_plot!(mplot, cfg, m.axis, data_processed, bounds, (tile, m.crs)) + end + + if haskey(m.current_tiles, tile) + move_in_front!(mplot, abs(m.zoom[] - tile.z), bounds) + else + move_to_back!(mplot, abs(m.zoom[] - tile.z), bounds) + end + + # Always move new plots to the front + mplot.visible = true + get_postprocess(cfg)(mplot) + m.plots[key] = (mplot, tile, bounds) + return +end + +############################ +#### Elevation Data plotting +#### + +struct ElevationData + elevation::AbstractMatrix{<: Number} + color::AbstractMatrix{<: Colorant} + elevation_range::Vec2d +end + +function Base.map(f::Function, data::ElevationData) + return ElevationData( + map(f, data.elevation), + data.color, + data.elevation_range + ) +end + +function get_bounds(tile::Tile, data::ElevationData, crs) + ext = MapTiles.extent(tile, crs) + mini, maxi = extrema(data.elevation) + origin = Vec3d(ext.X[1], ext.Y[1], mini) + w = Vec3d(ext.X[2] - ext.X[1], ext.Y[2] - ext.Y[1], maxi - mini) + return Rect3d(origin, w) +end + +function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::ElevationData, bounds::Rect, tile_crs) + # not so elegant with empty array, we may want to make this a bit nicer going forward + color = isempty(data.color) ? (;) : (color=data.color,) + mini, maxi = extrema(bounds) + p = Makie.surface!( + axis.scene, + (mini[1], maxi[1]), (mini[2], maxi[2]), data.elevation; + color..., + shading=Makie.NoShading, + inspectable=false, + colorrange=data.elevation_range, + config.attributes... + ) + return p +end + +function update_tile_plot!(plot::Surface, ::PlotConfig, ::AbstractAxis, data::ElevationData, bounds::Rect, tile_crs) + mini, maxi = extrema(bounds) + plot.args[1].val = (mini[1], maxi[1]) + plot.args[2].val = (mini[2], maxi[2]) + plot[3] = data.elevation + if !isempty(data.color) + plot.color = data.color + end + return +end + +############################ +#### Image Data plotting +#### + +const ImageData = AbstractMatrix{<:Colorant} + +function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) + mini, maxi = extrema(bounds) + plot = Makie.image!( + axis.scene, + (mini[1], maxi[1]), (mini[2], maxi[2]), rotr90(data); + inspectable=false, + config.attributes... + ) + return plot +end + +function update_tile_plot!(plot::Makie.Image, ::PlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) + mini, maxi = extrema(bounds) + plot[1] = (mini[1], maxi[1]) + plot[2] = (mini[2], maxi[2]) + plot[3] = rotr90(data) + return +end + + +############################ +#### PointCloudData Data plotting +#### + +struct PointCloudData + points::AbstractVector{<:Point3} + color::AbstractVector{<:Union{Colorant, Number}} + bounds::Rect3d + msize::Float32 +end + +get_bounds(::Tile, data::PointCloudData, crs) = data.bounds + +function Base.map(f::Function, data::PointCloudData) + return PointCloudData( + map(f, data.points), + data.color, + data.bounds, + data.msize + ) +end + +function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::PointCloudData, ::Rect, tile_crs) + p = Makie.scatter!( + axis.scene, data.points; + color=data.color, + marker=Makie.FastPixel(), + markersize=data.msize, + markerspace=:data, + fxaa=false, + inspectable=false, + config.attributes... + ) + return p +end + +function update_tile_plot!(plot::Makie.Scatter, ::PlotConfig, ::AbstractAxis, data::PointCloudData, bounds::Rect, tile_crs) + plot.color.val = data.color + plot[1] = data.points + plot.markersize = data.msize + return +end + +struct MeshScatterPlotconfig <: AbstractPlotConfig + plcfg::PlotConfig +end +MeshScatterPlotconfig(args...; attr...) = MeshScatterPlotconfig(PlotConfig(args...; attr...)) +get_preprocess(config::AbstractPlotConfig) = get_preprocess(config.plcfg) +get_postprocess(config::AbstractPlotConfig) = get_postprocess(config.plcfg) + +function create_tileplot!(config::MeshScatterPlotconfig, axis::AbstractAxis, data::PointCloudData, ::Rect, tile_crs) + m = Rect3f(Vec3f(0), Vec3f(1)) + p = Makie.meshscatter!( + axis.scene, data.points; + color=data.color, + marker=m, + markersize=data.msize, + inspectable=false, + config.plcfg.attributes... + ) + return p +end + +function update_tile_plot!(plot::Makie.MeshScatter, ::MeshScatterPlotconfig, ::AbstractAxis, data::PointCloudData, bounds::Rect, tile_crs) + plot.color = data.color + plot[1] = data.points + plot.markersize = data.msize + return +end + + + +############################ +#### Debug tile plotting (image only for now) +#### + +struct DebugPlotConfig <: AbstractPlotConfig + attributes::Dict{Symbol,Any} +end + +DebugPlotConfig(; plot_attributes...) = DebugPlotConfig(Dict{Symbol,Any}(plot_attributes)) + +function create_tileplot!(config::DebugPlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) + plot = Makie.poly!( + axis.scene, + bounds; + color=reverse(data; dims=1), + strokewidth=2, + strokecolor=:black, + inspectable=false, + stroke_depth_shift=-0.01f0, + config.attributes... + ) + return plot +end + +function update_tile_plot!(plot::Makie.Poly, ::DebugPlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) + plot[1] = bounds + plot.color = reverse(data; dims=1) + return +end diff --git a/src/tiles.jl b/src/tiles.jl new file mode 100644 index 0000000..cda9e0e --- /dev/null +++ b/src/tiles.jl @@ -0,0 +1,128 @@ +struct TileCache{TileFormat,Downloader} + provider::AbstractProvider + # Nothing for unavailable tiles, so we don't download them again + fetched_tiles::LRU{String,Union{Nothing, TileFormat}} + tile_queue::Channel{Tile} + # We also need to put! nothing for unavailable tiles into downloaded_tiles, + # so `Map` can clean up the expected tiles + downloaded_tiles::Channel{Tuple{Tile,Union{Nothing, TileFormat}}} + downloader::Vector{Downloader} +end + +get_tile_format(provider) = Matrix{RGB{N0f8}} + +get_downloader(provider) = ByteDownloader() + +function take_last!(c::Channel) + lock(c) + try + while isempty(c.data) + Base.check_channel_state(c) + wait(c.cond_take) + end + # function taken from Base.take_buffered, with just this line replaced to use `pop!` instead of `popfirst!` + v = pop!(c.data) + Base._increment_n_avail(c, -1) + notify(c.cond_put, nothing, false, false) # notify only one, since only one slot has become available for a put!. + return v + finally + unlock(c) + end +end + + +function run_loop(dl, tile_queue, fetched_tiles, provider, downloaded_tiles) + while isopen(tile_queue) || isready(tile_queue) + tile = take_last!(tile_queue) # priorize newly arrived tiles + result = nothing + try + @debug("downloading tile on thread $(Threads.threadid())") + # For providers which have to map the same data to different tiles + # Or providers that have e.g. additional paramers like date + # the url is a much better key than the tile itself + key = tile_key(provider, tile) + # if the provider knows it doesn't have a tile, it can return nothing + isnothing(key) && continue + result = get!(fetched_tiles, key) do + try + return fetch_tile(provider, dl, tile) + catch e + if isa(e, RequestError) + status = e.response.status + if (status == 404 || status == 500) + @warn "tile $(tile) not available, will not download again" maxlog = 10 + return nothing + end + end + rethrow(e) + end + end + catch e + @warn "Error while fetching tile on thread $(Threads.threadid())" exception = (e, catch_backtrace()) + nothing + end + put!(downloaded_tiles, (tile, result)) + yield() + end +end + +function Base.close(tiles::TileCache) + close(tiles.tile_queue) + close(tiles.downloaded_tiles) + empty!(tiles.fetched_tiles) +end + +function TileCache(provider; cache_size_gb=5, max_parallel_downloads=1) + TileFormat = get_tile_format(provider) + downloader = [get_downloader(provider) for i in 1:max_parallel_downloads] + fetched_tiles = LRU{String,Union{Nothing, TileFormat}}(; maxsize=cache_size_gb * 10^9, by=Base.summarysize) + downloaded_tiles = Channel{Tuple{Tile,Union{Nothing, TileFormat}}}(Inf) + tile_queue = Channel{Tile}(Inf) + async = Threads.nthreads(:default) <= 1 + if async && max_parallel_downloads > 1 + @warn "Multiple download threads are not supported with Threads.nthreads()==1, falling back to async. Start Julia with more threads for parallel downloads." + async = true + end + @assert max_parallel_downloads > 0 + for thread in 1:max_parallel_downloads + dl = downloader[thread] + if async + @async run_loop(dl, tile_queue, fetched_tiles, provider, downloaded_tiles) + else + Threads.@spawn run_loop(dl, tile_queue, fetched_tiles, provider, downloaded_tiles) + end + end + return TileCache{TileFormat,eltype(downloader)}(provider, fetched_tiles, tile_queue, downloaded_tiles, downloader) +end + +function fetch_tile(provider::AbstractProvider, downloader::AbstractDownloader, tile::Tile) + url = TileProviders.geturl(provider, tile.x, tile.y, tile.z) + isnothing(url) && return nothing + data = download_tile_data(downloader, provider, url) + return load_tile_data(provider, data) +end + +function load_tile_data(::AbstractProvider, downloaded::AbstractVector{UInt8}) + io = IOBuffer(downloaded) + format = FileIO.query(io) # this interrogates the magic bits to see what file format it is (JPEG, PNG, etc) + return FileIO.load(format) # this works because we have ImageIO loaded +end + +function Base.wait(tiles::TileCache; timeout=50) + # wait for all tiles to get downloaded + items = lock(tiles.tile_queue) do + copy(tiles.tile_queue.data) + end + tile_keys = filter!(!isnothing, map(t-> tile_key(tiles.provider, t), items)) + start = time() + while true + if isempty(tiles.tile_queue) && all(tk -> haskey(tiles.fetched_tiles, tk), tile_keys) + break + end + if time() - start > timeout + @warn "Timeout while waiting for tiles to download" + break + end + sleep(0.01) + end +end diff --git a/src/tyler-cam3d.jl b/src/tyler-cam3d.jl new file mode 100644 index 0000000..a11eade --- /dev/null +++ b/src/tyler-cam3d.jl @@ -0,0 +1,182 @@ +using Makie: ray_assisted_pick, @extract, mouseposition_px, parent_scene, screen_relative + +function translate_cam!(scene, cam::Camera3D, shift::Vec3d) + fb = cam.lookat[] .- cam.eyeposition[] + dir_left_right = normalize(cross(fb, cam.upvector[])) .* shift[1] + dir_front_back = normalize(Vec3(fb[1], fb[2], 0)) .* shift[2] + dir = (dir_front_back + dir_left_right) ./ 2 + cam.eyeposition[] += dir + cam.lookat[] += Vec3(dir[1], dir[2], 0) + update_cam!(scene, cam) +end + + +function _zoom!(scene, cam::Camera3D, zoom_step) + s = sign(1 - zoom_step) + eyepos = cam.eyeposition[] + lookat = cam.lookat[] + # just zoom up down + move_dir = Vec3(0, 0, 1) # -z + falloff = eyepos[3] / 10 + z_to_move = s .* zoom_step .* falloff + new_z = max(30, eyepos[3] - z_to_move) + cam.eyeposition[] = Vec3(eyepos[1], eyepos[2], new_z) + + # Calculate the direction vector from eyeposition to lookat + direction = Vec3(eyepos[1], eyepos[2], 0.0) .- lookat + ldistance = norm(direction) + lookat2d = lookat .+ ((ldistance / 15) .* s .* normalize(direction)) + cam.lookat[] = Vec3(lookat2d[1], lookat2d[2], 0) + + update_cam!(scene, cam) + return +end + +function signed_angle_between(v1::Vec3, v2::Vec3, normal::Vec3=Vec3d(0, 0, 1)) + # Normalize both vectors + v1n = normalize(v1) + v2n = normalize(v2) + + # Calculate the dot product + d = dot(v1n, v2n) + + # Calculate the cosine of the angle + cos_theta = clamp(d, -1.0, 1.0) + + # Calculate the angle in radians + angle_rad = acos(cos_theta) + + # Calculate the cross product + cross_product = cross(v1n, v2n) + if dot(cross_product, normal) < 0 + return -angle_rad + end + return angle_rad +end + + +function t_rotate_cam!(scene, cam::Camera3D, angles::VecTypes, from_mouse=false) + @extractvalue cam.controls (fix_x_key, fix_y_key, fix_z_key) + @extractvalue cam.settings (fixed_axis, circular_rotation, rotation_center) + + # This applies rotations around the x/y/z axis of the camera coordinate system + # x expands right, y expands up and z expands towards the screen + lookat = cam.lookat[] + eyepos = cam.eyeposition[] + up = cam.upvector[] # +y + viewdir = lookat - eyepos # -z + right = cross(viewdir, up) # +x + x_axis = right + y_axis = Vec3d(0, 0, 1) + z_axis = -viewdir + + rotation = Quaternionf(0, 0, 0, 1) + # restrict total quaternion rotation to one axis + rotation *= qrotation(y_axis, angles[2]) + rotation *= qrotation(x_axis, angles[1]) + rotation *= qrotation(z_axis, angles[3]) + new_viewdir = rotation * viewdir + new_angle = signed_angle_between(new_viewdir, Vec3d(0, 0, -1), -right) + if new_angle > 0 && new_angle < deg2rad(80) + new_lookat = ray_plane_intersection(Point3d(0), Vec3d(0, 0, 1), Point(eyepos), Vec(new_viewdir)) + isnothing(new_lookat) && return + cam.upvector[] = rotation * up + cam.lookat[] = new_lookat + update_cam!(scene, cam) + end + return +end + +function add_tyler_mouse_controls!(scene, cam::Camera3D) + @extract cam.controls (translation_button, rotation_button, reposition_button, scroll_mod) + @extract cam.settings ( + mouse_translationspeed, mouse_rotationspeed, mouse_zoomspeed, + cad, projectiontype, zoom_shift_lookat + ) + + last_mousepos = Base.RefValue(Vec2d(0, 0)) + dragging = Base.RefValue((false, false)) # rotation, translation + + e = events(scene) + + function compute_diff(delta) + if projectiontype[] == Makie.Perspective + # TODO wrong scaling? :( + ynorm = 2 * norm(cam.lookat[] - cam.eyeposition[]) * tand(0.5 * cam.fov[]) + return ynorm / size(scene, 2) * delta + else + viewnorm = norm(cam.eyeposition[] - cam.lookat[]) + return 2 * viewnorm / size(scene, 2) * delta + end + end + + # drag start/stop + on(camera(scene), e.mousebutton) do event + # Drag start translation/rotation + if event.action == Mouse.press && is_mouseinside(scene) + if ispressed(scene, translation_button[]) + last_mousepos[] = mouseposition_px(scene) + dragging[] = (false, true) + return Consume(true) + elseif ispressed(scene, rotation_button[]) + last_mousepos[] = mouseposition_px(scene) + dragging[] = (true, false) + return Consume(true) + end + # drag stop & repostion + elseif event.action == Mouse.release + consume = false + + # Drag stop translation/rotation + if dragging[][1] + mousepos = mouseposition_px(scene) + diff = compute_diff(last_mousepos[] .- mousepos) + last_mousepos[] = mousepos + dragging[] = (false, false) + translate_cam!(scene, cam, mouse_translationspeed[] .* Vec3d(diff[1], diff[2], 0.0)) + consume = true + elseif dragging[][2] + mousepos = mouseposition_px(scene) + dragging[] = (false, false) + rot_scaling = mouse_rotationspeed[] * (e.window_dpi[] * 0.005) + mp = (last_mousepos[] .- mousepos) .* 0.01 .* rot_scaling + last_mousepos[] = mousepos + t_rotate_cam!(scene, cam, Vec3d(-mp[2], mp[1], 0.0), true) + consume = true + end + + return Consume(consume) + end + + return Consume(false) + end + + # in drag + on(camera(scene), e.mouseposition) do mp + if dragging[][2] && ispressed(scene, translation_button[]) + mousepos = screen_relative(scene, mp) + diff = compute_diff(last_mousepos[] .- mousepos) + last_mousepos[] = mousepos + translate_cam!(scene, cam, mouse_translationspeed[] * Vec3d(diff[1], diff[2], 0.0)) + return Consume(true) + elseif dragging[][1] && ispressed(scene, rotation_button[]) + mousepos = screen_relative(scene, mp) + rot_scaling = mouse_rotationspeed[] * (e.window_dpi[] * 0.005) + mp = (last_mousepos[] .- mousepos) * 0.01 * rot_scaling + last_mousepos[] = mousepos + t_rotate_cam!(scene, cam, Vec3d(-mp[2], mp[1], 0.0), true) + return Consume(true) + end + return Consume(false) + end + + #zoom + on(camera(scene), e.scroll) do scroll + if is_mouseinside(scene) && ispressed(scene, scroll_mod[]) + zoom_step = (1.0 + 0.1 * mouse_zoomspeed[])^-scroll[2] + _zoom!(scene, cam, zoom_step) + return Consume(true) + end + return Consume(false) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 167be4b..58549b5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,26 +6,26 @@ using GeoInterface # Default london = Rect2f(-0.0921, 51.5, 0.04, 0.025) -m = wait(Tyler.Map(london; scale=1)) # waits until all tiles are displayed -@test isempty(m.tiles_being_added) -@test isempty(m.queued_but_not_downloaded) -@test length(m.displayed_tiles) == 24 -@test length(m.fetched_tiles) == 24 +m = Tyler.Map(london); m.figure.scene +s = display(m) # waits until all tiles are displayed +@test isempty(m.tiles.tile_queue) +@test length(m.current_tiles) == 25 +@test length(m.tiles.fetched_tiles) == 48 london = Rect2f(-0.0921, 51.5, 0.04, 0.025) -m = wait(Tyler.Map(london; scale=1, provider = Tyler.TileProviders.Google(), crs=Tyler.MapTiles.WGS84())) # waits until all tiles are displayed -@test isempty(m.tiles_being_added) -@test isempty(m.queued_but_not_downloaded) -@test length(m.displayed_tiles) == 24 -@test length(m.fetched_tiles) == 24 +m = Tyler.Map(london; scale=1, provider=Tyler.TileProviders.Google(), crs=Tyler.MapTiles.WGS84()) # waits until all tiles are displayed +s = display(m) # waits until all tiles are displayed +@test isempty(m.tiles.tile_queue) +@test length(m.current_tiles) == 35 +@test length(m.tiles.fetched_tiles) == 66 # test Extent input london = Extents.Extent(X=(-0.0921, -0.0521), Y=(51.5, 51.525)) -m = wait(Tyler.Map(london; scale=1)) # waits until all tiles are displayed -@test isempty(m.tiles_being_added) -@test isempty(m.queued_but_not_downloaded) -@test length(m.displayed_tiles) == 24 -@test length(m.fetched_tiles) == 24 +m = Tyler.Map(london; scale=1) # waits until all tiles are displayed +display(m) +@test isempty(m.tiles.tile_queue) +@test length(m.current_tiles) == 25 +@test length(m.tiles.fetched_tiles) == 48 @testset "Interfaces" begin from = Tyler.MapTiles.WebMercator() diff --git a/test/test-providers.jl b/test/test-providers.jl new file mode 100644 index 0000000..d29d032 --- /dev/null +++ b/test/test-providers.jl @@ -0,0 +1,194 @@ +using GeometryBasics, GLMakie, Tyler, TileProviders +using Tyler: ElevationProvider, GeoTilePointCloudProvider +using MapTiles +using MeshIO, FileIO +# https://api.3dbag.nl/api.html +# https://docs.3dbag.nl/en/delivery/webservices/ +amsti = load("amsterdam/amsterdam.obj") +tex = load("amsterdam/texture.png") +f, ax, pl = mesh(amsti; color=:gray) +boundingbox(pl) + +lat, lon = (52.395593, 4.884704) +delta = 0.01 +ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) +extrema(ext) +begin + lat, lon = (52.395593, 4.884704) + delta = 0.01 + ext = Rect2f(lon-delta/2, lat-delta/2, delta, delta) + m1 = Tyler.Map(ext; scale=0.5) +end +for (k, (pl, t, bb)) in m1.plots + linesegments!(m1.axis.scene, bb, color=:red, depth_shift=-0.1f0) +end + + +begin + lat, lon = (52.395593, 4.884704) + delta = 0.01 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + m1 = Tyler.Map(ext; max_parallel_downloads=1) + display(m1.figure.scene) +end + +begin + # lat, lon = (40.697211, -74.037523) + # lat, lon = (52.377428, 4.898387) + # lat, lon = (53.208252, 5.169539) + # lat, lon = (55.648466, 12.566546) + lat, lon = (47.087441, 13.377214) + delta = 0.3 + ext = Rect2f(lon-delta/2, lat-delta/2, delta, delta) + m = Tyler.Map3D(ext; provider=ElevationProvider()) +end + +begin + lat, lon = (47.087441, 13.377214) + delta = 0.3 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + cfg = Tyler.PlotConfig(preprocess=pc -> map(p -> p .* 2, pc), shading=FastShading, material=mat, colormap=:alpine) + m = Tyler.Map3D(ext; provider=ElevationProvider(nothing), plot_config=cfg) +end + +for i in 1:1000 + rotate_cam!(m.axis.scene, 0, 0.01, 0) + sleep(1/60) +end + +begin + lat, lon = (52.40459835, 4.84763329) + delta = 0.006 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + subset="AHN1_T" # Takes reasonably long to load (~1-5mb compressed per tile) + # subset="AHN4_T" # Takes _really_ long to load, even from disk (~300mb compressed points per tile) + provider = GeoTilePointCloudProvider(subset=subset) + image = ElevationProvider(nothing) + cfg = Tyler.MeshScatterPlotconfig(markersize=5, marker=Rect3f(Vec3f(0), Vec3f(1))) + m1 = Tyler.Map3D(ext; provider=provider, plot_config=cfg, max_parallel_downloads=1) + cfg = Tyler.PlotConfig(preprocess=pc -> map(p -> p .* 2, pc), shading=FastShading, colormap=:alpine, postprocess=(p-> translate!(p, 0, 0, -1f0))) + m2 = Tyler.Map3D(m1; provider=image, plot_config=cfg, max_parallel_downloads=1) + m1 +end + +using GeometryBasics, GLMakie, Tyler, TileProviders +using Tyler: ElevationProvider, GeoTilePointCloudProvider +using MapTiles +using MeshIO, FileIO +begin + lat, lon = (52.395593, 4.884704) + delta = 0.02 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + subset = "AHN1_T" # Takes reasonably long to load (~1-5mb compressed per tile) + # subset = "AHN4_T" # Takes _really_ long to load, even from disk (~300mb compressed points per tile) + provider = GeoTilePointCloudProvider(subset=subset) + m1 = Tyler.Map3D(ext; provider=provider, max_parallel_downloads=1) + cfg = Tyler.PlotConfig(shading=FastShading, colormap=:alpine, postprocess=(p -> translate!(p, 0, 0, -100.0f0))) + m2 = Tyler.Map3D(ext; provider=ElevationProvider(), figure=m1.figure, axis=m1.axis, max_parallel_downloads=1, plot_config=cfg) + m1 +end +max_zoom(m1) +Tyler.approx_tiles(m1, , 1000) + +m1.current_tiles +m1.tiles.tile_queue +m1.plots +m1.should_get_plotted + +Tyler.cleanup_queue!(m1, Tyler.OrderedSet{Tile}()) + +m1.axis.scene.camera_controls.near[] = 100 +m1.axis.scene.camera_controls.far[] = 10000 +update_cam!(m1.axis.scene) +pp1 = Makie.project.((m1.axis.scene,), :data, :clip, maximum.(bb1)) +pp2 = Makie.project.((m1.axis.scene,), :data, :clip, maximum.(bb2)) +mean(Float32.(abs.(last.(pp1) .- last.(pp2[1:12])))) + +begin + lat, lon = (52.395593, 4.884704) + delta = 0.02 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + subset = "AHN1_T" # Takes reasonably long to load (~1-5mb compressed per tile) + # subset = "AHN4_T" # Takes _really_ long to load, even from disk (~300mb compressed points per tile) + provider = GeoTilePointCloudProvider(subset=subset) + m1 = Tyler.Map3D(ext; provider=provider) + m2 = Tyler.Map3D(ext; provider=ElevationProvider(), figure=m1.figure, axis=m1.axis) + m1 +end + +using RPRMakie, FileIO + +function plastic_material() + return (type=:Uber, reflection_color=Vec4f(1), + reflection_weight=Vec4f(1), reflection_roughness=Vec4f(0.1), + reflection_anisotropy=Vec4f(0), reflection_anisotropy_rotation=Vec4f(0), + reflection_metalness=Vec4f(0), reflection_ior=Vec4f(1.4), refraction_weight=Vec4f(0), + coating_weight=Vec4f(0), sheen_weight=Vec4f(0), emission_weight=Vec3f(0), + transparency=Vec4f(0), reflection_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), + emission_mode=UInt(RPR.RPR_UBER_MATERIAL_EMISSION_MODE_SINGLESIDED), + coating_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), sss_multiscatter=true, + refraction_thin_surface=true) +end +begin + lat, lon = (52.40459835, 4.84763329) + delta = 0.005 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + provider = GeoTilePointCloudProvider() + mat = plastic_material() + cfg = Tyler.MeshScatterPlotconfig(material=mat) + m = Tyler.Map3D(ext; figure=f, axis=ax, provider=provider, max_plots=3, plot_config=cfg) +end + + +using RPRMakie, FileIO +function render_rpr(m, name, radiance=1000000) + wait(m) + ax = m.axis + cam = ax.scene.camera_controls + lightpos = Vec3f(cam.lookat[][1], cam.lookat[][2], cam.eyeposition[][3]) + lights = [ + EnvironmentLight(1.5, load(RPR.assetpath("studio026.exr"))), + PointLight(lightpos, RGBf(radiance, radiance * 0.9, radiance * 0.9)) + ] + empty!(ax.scene.lights) + append!(ax.scene.lights, lights) + save("$(name).png", ax.scene; plugin=RPR.Northstar, backend=RPRMakie, iterations=2000) +end + +function plastic_material() + return (type=:Uber, reflection_color=Vec4f(1), + reflection_weight=Vec4f(1), reflection_roughness=Vec4f(0.1), + reflection_anisotropy=Vec4f(0), reflection_anisotropy_rotation=Vec4f(0), + reflection_metalness=Vec4f(0), reflection_ior=Vec4f(1.4), refraction_weight=Vec4f(0), + coating_weight=Vec4f(0), sheen_weight=Vec4f(0), emission_weight=Vec3f(0), + transparency=Vec4f(0), reflection_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), + emission_mode=UInt(RPR.RPR_UBER_MATERIAL_EMISSION_MODE_SINGLESIDED), + coating_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), sss_multiscatter=true, + refraction_thin_surface=true) +end + +begin + lat, lon = (47.087441, 13.377214) + delta = 0.5 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + mat = (type=:Uber, roughness=0.2, ior=1.390) + cfg = Tyler.PlotConfig(preprocess=pc -> map(p -> p .* 2, pc), shading=FastShading, material=plastic, colormap=:alpine) + m = Tyler.Map3D(ext; provider=ElevationProvider(nothing), plot_config=cfg, max_plots=5) + render_rpr(m, "alpine", 10000000) +end + +begin + lat, lon = (52.40459835229174, 4.84763329882317) + delta = 0.005 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + subset = "AHN1_T" # Takes reasonably long to load (~1-5mb compressed per tile) + # subset = "AHN4_T" # Takes _really_ long to load, even from disk (~300mb compressed points per tile) + provider = GeoTilePointCloudProvider(subset=subset) + mat = plastic_material() + cfg = Tyler.MeshScatterPlotconfig(markersize=5, material=mat) + m = Tyler.Map3D(ext; provider=provider, plot_config=cfg, max_plots=3, size=(2000, 2000)) + cfg = Tyler.PlotConfig(preprocess=pc -> map(p -> p .* 2, pc), shading=FastShading, material=mat, colormap=:Blues) + m2 = Tyler.Map3D(m; provider=ElevationProvider(nothing), plot_config=cfg, max_plots=5) + wait(m) + render_rpr(m, "pointclouds") +end diff --git a/test/tests.jl b/test/tests.jl new file mode 100644 index 0000000..4f79356 --- /dev/null +++ b/test/tests.jl @@ -0,0 +1,12 @@ +begin + lat, lon = (52.395593, 4.884704) + delta = 0.02 + ext = Rect2f(lon - delta / 2, lat - delta / 2, delta, delta) + subset = "AHN1_T" # Takes reasonably long to load (~1-5mb compressed per tile) + # subset = "AHN4_T" # Takes _really_ long to load, even from disk (~300mb compressed points per tile) + provider = GeoTilePointCloudProvider(subset=subset) + m1 = Tyler.Map3D(ext; provider=provider) + wait(m1) + unique_plots = unique(Tyler.tile_key.((m1.provider,), keys(m1.current_tiles))) + @test length(unique_plots) == length(m1.plots) +end