diff --git a/src/Tyler.jl b/src/Tyler.jl index f905ab67..9fd8c6a3 100644 --- a/src/Tyler.jl +++ b/src/Tyler.jl @@ -1,30 +1,33 @@ module Tyler +using ArchGDAL +using Downloads +using FileIO +using GeometryBasics +using ImageIO +using LinearAlgebra +using Makie +using PointClouds +using Proj +using Scratch +using Statistics + +import GeoFormatTypes as GFT +import GeometryOps as GO + 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, Vec3d, Point2d, Point3d, Point4d +using GeometryBasics: GeometryBasics, GLTriangleFace, Point2f, Vec2f, Rect2f, Rect2, Rect, Vec3d, Point2d, Point3d, Point4d +using GeometryBasics: decompose, decompose_uv using HTTP: HTTP using LRUCache: LRUCache, LRU +using Makie: AbstractAxis, Observable, Figure, Axis, LScene, RGBAf +using Makie: on, isopen, meta, mesh!, translate!, scale!, Plot using MapTiles: MapTiles, Tile, TileGrid, web_mercator, wgs84, CoordinateReferenceSystemFormat -using Makie: Makie, Observable, Figure, Axis, LScene, RGBAf, on, isopen, 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 -using PointClouds -using ArchGDAL -import GeoFormatTypes as GFT -using Downloads -using Scratch -using ImageIO, FileIO -import GeometryOps as GO - const CACHE_PATH = Ref("") @@ -38,15 +41,15 @@ abstract type FetchingScheme end abstract type AbstractMap end include("downloader.jl") -include("tiles.jl") +include("tile-cache.jl") include("map.jl") -include("3d-map.jl") -include("tyler-cam3d.jl") +include("map3d.jl") +include("cam3d.jl") include("tile-plotting.jl") include("tile-fetching.jl") +include("provider/shared.jl") include("provider/interpolations.jl") include("provider/elevation/elevation-provider.jl") include("provider/pointclouds/geotiles-pointcloud-provider.jl") - end diff --git a/src/tyler-cam3d.jl b/src/cam3d.jl similarity index 99% rename from src/tyler-cam3d.jl rename to src/cam3d.jl index a11eadee..b9ee1d6b 100644 --- a/src/tyler-cam3d.jl +++ b/src/cam3d.jl @@ -179,4 +179,4 @@ function add_tyler_mouse_controls!(scene, cam::Camera3D) end return Consume(false) end -end +end \ No newline at end of file diff --git a/src/downloader.jl b/src/downloader.jl index bcea7c5b..de3efc38 100644 --- a/src/downloader.jl +++ b/src/downloader.jl @@ -9,7 +9,6 @@ struct ByteDownloader <: AbstractDownloader 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) @@ -31,7 +30,6 @@ struct PathDownloader <: AbstractDownloader 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) @@ -40,15 +38,13 @@ function PathDownloader(cache_dir; timeout=5, cache_size_gb=5) 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) + 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 + +_unique_filename(url) = string(hash(url)) diff --git a/src/map.jl b/src/map.jl index eccf2bca..822c842a 100644 --- a/src/map.jl +++ b/src/map.jl @@ -1,9 +1,3 @@ - -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 @@ -57,81 +51,6 @@ struct Map{Ax<:Makie.AbstractAxis} <: AbstractMap max_plots::Int scale::Float64 end - -function setup_figure_and_axis!(figure::Makie.Figure, axis, ext_target, crs) - setup_axis!(axis, ext_target, crs) - return figure, axis -end - -function setup_figure_and_axis!(figure::Makie.Figure, axis_kws_nt::NamedTuple, ext_target, crs) - axis_kws = Dict(pairs(axis_kws_nt)) - AxisType = pop!(axis_kws, :type, Axis) - - axis = AxisType(figure[1, 1]; axis_kws...) - - setup_axis!(axis, ext_target, crs) - - return figure, axis -end - -_get_parent_layout(gp::Makie.GridPosition) = _get_parent_layout(gp.layout) -_get_parent_layout(gp::Makie.GridSubposition) = _get_parent_layout(gp.layout) -_get_parent_layout(gl::Makie.GridLayout) = gl - - -_get_parent_figure(fig::Makie.Figure) = fig -_get_parent_figure(gl::Makie.GridLayout) = _get_parent_figure(gl.parent) -_get_parent_figure(gp::Makie.GridPosition) = _get_parent_figure(_get_parent_layout(gp.layout)) -_get_parent_figure(gp::Makie.GridSubposition) = _get_parent_figure(_get_parent_layout(gp.layout)) - -function setup_figure_and_axis!(figure::GridPosition, axis, ext_target, crs) - error(""" - You have tried to construct a `Map` at a given grid position, - but with a materialized axis of type $(typeof(axis)). - - You can only do this if you let Tyler construct the axis, - by passing its parameters as a NamedTuple - (like `axis = (; type = Axis, ...)`). - """) -end - -function setup_figure_and_axis!(gridposition::GridPosition, axis_kws_nt::NamedTuple, ext_target, crs) - figure = _get_parent_figure(gridposition) - - axis_kws = Dict(pairs(axis_kws_nt)) - AxisType = pop!(axis_kws, :type, Axis) - - axis = AxisType(gridposition; axis_kws...) - - setup_axis!(axis, ext_target, crs) - - return figure, axis -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. @@ -158,18 +77,6 @@ function Map(m::Map; kw...) 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.foreground_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), @@ -247,29 +154,68 @@ function Map(extent, extent_crs=wgs84; 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.foreground_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 +function setup_figure_and_axis!(figure::Makie.Figure, axis, ext_target, crs) + setup_axis!(axis, ext_target, crs) + return figure, axis +end +function setup_figure_and_axis!(figure::Makie.Figure, axis_kws_nt::NamedTuple, ext_target, crs) + axis_kws = Dict(pairs(axis_kws_nt)) + AxisType = pop!(axis_kws, :type, Axis) + + axis = AxisType(figure[1, 1]; axis_kws...) + + setup_axis!(axis, ext_target, crs) + + return figure, axis +end +function setup_figure_and_axis!(figure::GridPosition, axis, ext_target, crs) + error(""" + You have tried to construct a `Map` at a given grid position, + but with a materialized axis of type $(typeof(axis)). + + You can only do this if you let Tyler construct the axis, + by passing its parameters as a NamedTuple + (like `axis = (; type = Axis, ...)`). + """) +end +function setup_figure_and_axis!(gridposition::GridPosition, axis_kws_nt::NamedTuple, ext_target, crs) + figure = _get_parent_figure(gridposition) + + axis_kws = Dict(pairs(axis_kws_nt)) + AxisType = pop!(axis_kws, :type, Axis) + + axis = AxisType(gridposition; axis_kws...) + + setup_axis!(axis, ext_target, crs) + + return figure, axis +end + +_get_parent_layout(gp::Makie.GridPosition) = _get_parent_layout(gp.layout) +_get_parent_layout(gp::Makie.GridSubposition) = _get_parent_layout(gp.layout) +_get_parent_layout(gl::Makie.GridLayout) = gl + +_get_parent_figure(fig::Makie.Figure) = fig +_get_parent_figure(gl::Makie.GridLayout) = _get_parent_figure(gl.parent) +_get_parent_figure(gp::Makie.GridPosition) = _get_parent_figure(_get_parent_layout(gp.layout)) +_get_parent_figure(gp::Makie.GridSubposition) = _get_parent_figure(_get_parent_layout(gp.layout)) + +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 +toggle_visibility!(m::Map) = m.axis.scene.visible[] = !m.axis.scene.visible[] + function tile_reloader(map::Map{Axis}) axis = map.axis throttled = Makie.Observables.throttle(0.2, axis.finallimits) @@ -283,12 +229,29 @@ function plotted_tiles(m::Map) return getindex.(values(m.plots), 2) 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 + +# Package interface methods + 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 methods + Base.showable(::MIME"image/png", ::AbstractMap) = true function Base.show(io::IO, m::MIME"image/png", map::AbstractMap) @@ -301,17 +264,33 @@ function Base.display(map::AbstractMap) Base.display(map.figure.scene) end - -function remove_unused!(m::AbstractMap, tile::Tile) - return remove_unused!(m, tile_key(m.provider, tile)) +function Base.close(m::Map) + cleanup_queue!(m, OrderedSet{Tile}()) + empty!(m.foreground_tiles) + empty!(m.unused_plots) + empty!(m.plots) + empty!(m.should_get_plotted) + close(m.tiles) 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 +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 - return nothing -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.foreground_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 \ No newline at end of file diff --git a/src/3d-map.jl b/src/map3d.jl similarity index 91% rename from src/3d-map.jl rename to src/map3d.jl index 18d15092..171a96ee 100644 --- a/src/3d-map.jl +++ b/src/map3d.jl @@ -1,5 +1,29 @@ +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 +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 -using GeometryBasics, LinearAlgebra +# Tyler methods for LScene function setup_axis!(axis::LScene, ext_target, crs) # Disconnect all events @@ -21,24 +45,6 @@ function setup_axis!(axis::LScene, ext_target, crs) 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 @@ -50,6 +56,7 @@ function tile_reloader(m::Map{LScene}) end end +# 3d utils function frustum_snapshot(cam::Camera) bottom_left = Point4d(-1, -1, 1, 1) @@ -110,6 +117,7 @@ function frustrum_plane_intersection(cam::Camera, camc::Camera3D) end end +# TODO: move to GeometryBasics 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 diff --git a/src/provider/elevation/elevation-provider.jl b/src/provider/elevation/elevation-provider.jl index dc79c76a..9f0b4139 100644 --- a/src/provider/elevation/elevation-provider.jl +++ b/src/provider/elevation/elevation-provider.jl @@ -12,12 +12,6 @@ struct ElevationProvider <: 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) @@ -25,17 +19,22 @@ function ElevationProvider(provider=TileProviders.Esri(:WorldImagery); cache_siz ElevationProvider(provider, fetched_tiles, downloader) end +# TileProviders interface +TileProviders.options(::ElevationProvider) = nothing +TileProviders.min_zoom(::ElevationProvider) = 0 +TileProviders.max_zoom(::ElevationProvider) = 16 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 +# Tyler interface +get_tile_format(::ElevationProvider) = ElevationData +file_ending(::ElevationProvider) = ".lerc" 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"]) diff --git a/src/provider/interpolations.jl b/src/provider/interpolations.jl index 6d07c47d..2a24c380 100644 --- a/src/provider/interpolations.jl +++ b/src/provider/interpolations.jl @@ -14,14 +14,15 @@ 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) +# TileProviders interface function TileProviders.geturl(::Interpolator, x::Integer, y::Integer, z::Integer) return "$x,$y,$z" end +# Tyler interface get_downloader(::Interpolator) = NoDownload() function fetch_tile(interpolator::Interpolator, ::NoDownload, tile::Tile) @@ -30,6 +31,7 @@ function fetch_tile(interpolator::Interpolator, ::NoDownload, tile::Tile) return [_col(interpolator, i) for i in z] end + # 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)) @@ -48,4 +50,4 @@ function _tile2positions(x, y, z) lons = [_tile2lng(x + i, z) for i in rng, j in rng] lats = [_tile2lat(y + j, z) for i in rng, j in rng] return (lons, lats) -end +end \ No newline at end of file diff --git a/src/provider/pointclouds/geotiles-pointcloud-provider.jl b/src/provider/pointclouds/geotiles-pointcloud-provider.jl index 4dc2d781..b46911e9 100644 --- a/src/provider/pointclouds/geotiles-pointcloud-provider.jl +++ b/src/provider/pointclouds/geotiles-pointcloud-provider.jl @@ -1,4 +1,3 @@ - """ GeoTilePointCloudProvider(subset="AHN1_T") @@ -15,48 +14,15 @@ struct GeoTilePointCloudProvider <: TileProviders.AbstractProvider 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" - +# TileProviders interface +TileProviders.options(::GeoTilePointCloudProvider) = nothing +TileProviders.min_zoom(::GeoTilePointCloudProvider) = 16 +TileProviders.max_zoom(::GeoTilePointCloudProvider) = 16 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") @@ -64,11 +30,12 @@ function TileProviders.geturl(p::GeoTilePointCloudProvider, x::Integer, y::Integ 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 +# Tyler interface +get_tile_format(::GeoTilePointCloudProvider) = PointCloudData +file_ending(::GeoTilePointCloudProvider) = ".laz" +function get_downloader(::GeoTilePointCloudProvider) + cache_dir = joinpath(CACHE_PATH[], "GeoTilePointCloudProvider") + return PathDownloader(cache_dir) end function load_tile_data(provider::GeoTilePointCloudProvider, path::String) @@ -94,3 +61,36 @@ function load_tile_data(provider::GeoTilePointCloudProvider, path::String) bb = Rect3d(extrema[1], extrema[2] .- extrema[1]) return PointCloudData(points, color, bb, best_markersize[provider.subset]) end + +# Geotile utils + +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 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 \ No newline at end of file diff --git a/src/provider/shared.jl b/src/provider/shared.jl new file mode 100644 index 00000000..4e371a40 --- /dev/null +++ b/src/provider/shared.jl @@ -0,0 +1,22 @@ + +# Tyler provider defauls + +get_tile_format(provider) = Matrix{RGB{N0f8}} +get_downloader(provider) = ByteDownloader() + +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 tile_key(provider::AbstractProvider, tile::Tile) + return TileProviders.geturl(provider, tile.x, tile.y, tile.z) +end \ No newline at end of file diff --git a/src/tiles.jl b/src/tile-cache.jl similarity index 86% rename from src/tiles.jl rename to src/tile-cache.jl index 39cd17fc..b34f485e 100644 --- a/src/tiles.jl +++ b/src/tile-cache.jl @@ -8,10 +8,56 @@ struct TileCache{TileFormat,Downloader} downloaded_tiles::Channel{Tuple{Tile,Union{Nothing, TileFormat}}} downloader::Vector{Downloader} 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) -get_tile_format(provider) = Matrix{RGB{N0f8}} + 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 -get_downloader(provider) = ByteDownloader() +# Base methods + +function Base.close(tiles::TileCache) + close(tiles.tile_queue) + close(tiles.downloaded_tiles) + empty!(tiles.fetched_tiles) +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 function take_last!(c::Channel) lock(c) @@ -30,7 +76,6 @@ function take_last!(c::Channel) 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 @@ -64,66 +109,4 @@ function run_loop(dl, tile_queue, fetched_tiles, provider, downloaded_tiles) 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 +end \ No newline at end of file diff --git a/src/tile-fetching.jl b/src/tile-fetching.jl index d497c6ca..6fcf294c 100644 --- a/src/tile-fetching.jl +++ b/src/tile-fetching.jl @@ -1,3 +1,4 @@ +# Queue management function queue_plot!(m::Map, tile) key = tile_key(m.provider, tile) @@ -187,10 +188,18 @@ function get_tiles_for_area(m::Map{LScene}, ::Tiling3D, (cam, camc)::Tuple{Camer camc.far[] = maxdist camc.near[] = eyepos[3] * 0.01 update_cam!(m.axis.scene) - foreground = tiles_from_poly(m, points) + foreground = tiles_from_poly(m, points; zshift=0) background = OrderedSet{Tile}() + # for i in 2:2:6 + # tiles = tiles_from_poly(m, points; zshift=-i) + # union!(foreground, tiles) + # end offscreen = OrderedSet{Tile}() - return (; foreground, background, offscreen) + tiles = (; foreground, background, offscreen) + return tiles + # Reverse the order of the groups. Reversing the ranges + # above doesn't have the same effect due to then unions + # return map(OrderedSet ∘ reverse ∘ collect, tiles) end function get_tiles_for_area(m::Map{LScene}, s::SimpleTiling, (cam, camc)::Tuple{Camera,Camera3D}) area = area_around_lookat(camc) @@ -206,6 +215,7 @@ function get_resolution(map::Map) return isnothing(screen) ? size(map.axis.scene) .* 1.5 : size(screen) end +# TODO this will be in Extents.jl soon, so remove function grow_extent(area::Union{Rect,Extent}, factor) Extent(map(Extents.bounds(area)) do axis_bounds span = axis_bounds[2] - axis_bounds[1] @@ -225,7 +235,6 @@ function optimal_zoom(m::Map, diagonal) 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)) diff --git a/src/tile-plotting.jl b/src/tile-plotting.jl index 9ebc7683..b40c4419 100644 --- a/src/tile-plotting.jl +++ b/src/tile-plotting.jl @@ -47,7 +47,8 @@ extent = Extent(; X=(lon - delta / 2, lon + delta / 2), Y=(lat - delta / 2, lat 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) +PlotConfig(; preprocess=identity, postprocess=identity, plot_attributes...) = + PlotConfig(Dict{Symbol,Any}(plot_attributes), preprocess, postprocess) function get_bounds(tile::Tile, data, crs) @@ -193,7 +194,9 @@ function get_bounds(tile::Tile, data::ElevationData, crs) return Rect3d(origin, w) end -function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::ElevationData, bounds::Rect, tile_crs) +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) @@ -211,7 +214,9 @@ function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::Elevatio return p end -function update_tile_plot!(plot::Surface, ::PlotConfig, ::AbstractAxis, data::ElevationData, bounds::Rect, tile_crs) +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]) @@ -228,7 +233,9 @@ end const ImageData = AbstractMatrix{<:Colorant} -function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) +function create_tileplot!( + config::PlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs +) mini, maxi = extrema(bounds) plot = Makie.image!( axis.scene, @@ -240,7 +247,9 @@ function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::ImageDat return plot end -function update_tile_plot!(plot::Makie.Image, ::PlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) +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]) @@ -270,7 +279,9 @@ function Base.map(f::Function, data::PointCloudData) ) end -function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::PointCloudData, ::Rect, tile_crs) +function create_tileplot!( + config::PlotConfig, axis::AbstractAxis, data::PointCloudData, ::Rect, tile_crs +) p = Makie.scatter!( axis.scene, data.points; color=data.color, @@ -284,7 +295,9 @@ function create_tileplot!(config::PlotConfig, axis::AbstractAxis, data::PointClo return p end -function update_tile_plot!(plot::Makie.Scatter, ::PlotConfig, ::AbstractAxis, data::PointCloudData, bounds::Rect, tile_crs) +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 @@ -298,7 +311,9 @@ MeshScatterPlotconfig(args...; attr...) = MeshScatterPlotconfig(PlotConfig(args. 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) +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; @@ -311,7 +326,9 @@ function create_tileplot!(config::MeshScatterPlotconfig, axis::AbstractAxis, dat return p end -function update_tile_plot!(plot::Makie.MeshScatter, ::MeshScatterPlotconfig, ::AbstractAxis, data::PointCloudData, bounds::Rect, tile_crs) +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 @@ -326,10 +343,12 @@ end struct DebugPlotConfig <: AbstractPlotConfig attributes::Dict{Symbol,Any} end +DebugPlotConfig(; plot_attributes...) = + DebugPlotConfig(Dict{Symbol,Any}(plot_attributes)) -DebugPlotConfig(; plot_attributes...) = DebugPlotConfig(Dict{Symbol,Any}(plot_attributes)) - -function create_tileplot!(config::DebugPlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) +function create_tileplot!( + config::DebugPlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs +) plot = Makie.poly!( axis.scene, bounds; @@ -343,8 +362,10 @@ function create_tileplot!(config::DebugPlotConfig, axis::AbstractAxis, data::Ima return plot end -function update_tile_plot!(plot::Makie.Poly, ::DebugPlotConfig, axis::AbstractAxis, data::ImageData, bounds::Rect, tile_crs) +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 +end \ No newline at end of file diff --git a/test/test-providers.jl b/test/test-providers.jl index bc5fefdd..631ed6c7 100644 --- a/test/test-providers.jl +++ b/test/test-providers.jl @@ -110,7 +110,7 @@ begin # 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) + # m2 = Tyler.Map3D(ext; provider=ElevationProvider(), figure=m1.figure, axis=m1.axis) m1 end