View sattelite data in Makie.jl: Globe at lower zoom levels, only load data required for the current camera view

makie-globe-zoom.jl

using GeometryBasics, LinearAlgebra, GLMakie, FileIO, CoordinateTransformations, ColorSchemes

function intersect_line_sphere(p::Point3f, d::Point3f, c::Point3f, sphere_radius::Float32)
    # Use indexing to access point/vector components
    # Calculate coefficients of the quadratic equation (a*t^2 + b*t + c = 0)
    a = dot(d, d)
    b = 2 * dot(d, p .- c)
    c = dot(p .- c, p .- c) .- sphere_radius^2

    # Calculate discriminant
    discriminant = b^2 - 4 * a * c

    if discriminant < 0
        return nothing
    elseif discriminant == 0
        # One intersection point (tangent)
        t = -b / (2 * a)
        intersection_point = p + t * d
        return Point3f(intersection_point...)
    else
        # Two intersection points
        t1 = (-b + sqrt(discriminant)) / (2 * a)
        t2 = (-b - sqrt(discriminant)) / (2 * a)
        intersection_point1 = p + t1 * d
        intersection_point2 = p + t2 * d
        return Point3f(intersection_point1...), Point3f(intersection_point2...)
    end
end

function get_frustrum(cam)
    rect_ps = [
        Point3f(1, 1, 1), Point3f(-1, 1, 1),
        Point3f(1, -1, 1), Point3f(-1, -1, 1)]

    inv_pv = inv(cam.projectionview[])
    return map(rect_ps) do p
        p = inv_pv * to_ndim(Point4f, p, 1)
        return p[Vec(1, 2, 3)] / p[4]
    end
end

fig = Figure()
globe_scene = LScene(fig[1, 1]; show_axis=false)
globe_scene.scene.clear = true
globe_cam = Makie.Camera3D(
    globe_scene.scene,
    show_axis=false,
    projectiontype=Makie.Perspective,
    clipping_mode=:static,

    # disable translation to prevent skewing
    mouse_translationspeed=0,
    keyboard_translationspeed=0
)
globe_cam.near[] = 0.0001f0

earth_texture = load(Makie.assetpath("earth.png"))
earth = mesh!(globe_scene, Sphere(Point3f(0), 1.0f0), color=earth_texture)

center!(globe_scene.scene)

cam = globe_scene.scene.camera
eyeposition = globe_cam.eyeposition
lookat = globe_cam.lookat
frustrum = map(pv -> get_frustrum(cam), cam.projectionview)

cam_scene = LScene(fig[1, 2])
center!(cam_scene.scene)
meshscatter!(cam_scene, frustrum, color=:blue)

intersections = map(frustrum, eyeposition) do frustrum, eypos
    result = Point3f[]
    for point in frustrum
        res = intersect_line_sphere(Point(eypos), Point(eypos .- point), Point3f(0), 1.0f0)
        if !isnothing(res)
            if res isa Tuple
                res = sort!([res...], by=(x -> norm(x - eypos)))[1]
            end
            push!(result, res)
        end
    end
    return result
end

visible = map(intersections) do intersections
    if length(intersections) == 4
        earth.visible = false
    else
        earth.visible = true
    end
end

plane = map(intersections) do points
    if length(points) == 4
        return GeometryBasics.Mesh(meta(points;
                normals=fill(Vec3f(0), 4),
                uv=[Vec2f(0, 1), Vec2f(1, 1), Vec2f(0, 0), Vec2f(1, 0)]
            ),
            [GLTriangleFace(1, 2, 3), GLTriangleFace(2, 3, 4)])
    else
        return GeometryBasics.Mesh(meta([Point3f(0)]; normals=[Vec3f(0)], uv=[Vec2f(NaN)]), GLTriangleFace[])
    end
end

sphere_to_cart_trans = SphericalFromCartesian()
plane_bbox = map(intersections) do points
    if length(points) == 4
        # move date line from meridian to common date line
        p = sphere_to_cart_trans.(points) .|> x -> (x.θ > 0 ? x.θ - π : x.θ + π, x.ϕ) .|> rad2deg
        return p
    else
        return [(-180.0, 180.0), (180.0, -180.0), (-180.0, -180.0), (180.0, 180.0)]
    end
end

# TODO: Replace mock texture with actual data inside the calculated bounding box
s = 100
values = reshape(range(1, 255, s * s), (s, s)) .|> round .|> Integer
plane_texture = map(x -> ColorSchemes.viridis[x], values)

scatter!(cam_scene, eyeposition, color=:black)
meshscatter!(cam_scene, intersections, color=:white)
wireframe!(cam_scene, Sphere(Point3f(0), 1.0f0))
mesh!(cam_scene, plane, color=plane_texture, shading=NoShading)
mesh!(globe_scene, plane, color=plane_texture, shading=NoShading)
scatter(fig[1, 3], plane_bbox; axis=(ylabel="latitude [°]", xlabel="longitude [°]", title="Bounding Box"))

fig

TODO: Transform a view (current frustrum) of the global data array to the desired bounding box. E.g. via skewing/shearing?

The frustum forms a plane that intersects with the sphere resulting into a circle map. However, the viewport is rectangular. We can fill corner points with those being even further away e.g. using an orthographic projection.

• The corner points of the frustrum form quads that aren't even trapezoids. The texture pixels would be very irregular.
• For instance, Google Earth uses a perspective projection i.e, the virtual camera looking at the globe. Not very useful to produce equal area, angle or distance map.