From 54d147329183c67fd8ca30d0c020f9db4c8a6cfc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 29 Oct 2024 17:40:28 -0700 Subject: [PATCH 01/12] Amend examples + add globeaxis framework --- examples/earth3d.jl | 30 +++ examples/geodesy_raster_animation.jl | 52 ++++ examples/piracy_at_sea.jl | 100 +++++++ examples/specialized/true_icosphere.jl | 253 ++++++++++++++++++ examples/usa_inset.jl | 94 +++++++ src/sphere/globeaxis.jl | 356 +++++++++++++++++++++++++ src/sphere/icosphere.jl | 210 +++++++++++++++ src/sphere/unit_sphere_transforms.jl | 55 ++++ 8 files changed, 1150 insertions(+) create mode 100644 examples/earth3d.jl create mode 100644 examples/geodesy_raster_animation.jl create mode 100644 examples/piracy_at_sea.jl create mode 100644 examples/specialized/true_icosphere.jl create mode 100644 examples/usa_inset.jl create mode 100644 src/sphere/globeaxis.jl create mode 100644 src/sphere/icosphere.jl create mode 100644 src/sphere/unit_sphere_transforms.jl diff --git a/examples/earth3d.jl b/examples/earth3d.jl new file mode 100644 index 00000000..6d2f844c --- /dev/null +++ b/examples/earth3d.jl @@ -0,0 +1,30 @@ +using GLMakie, GeoMakie, Geodesy + +transf = Geodesy.ECEFfromLLA(Geodesy.WGS84()) + +f, a, p = meshimage(-180..180, -90..90, GeoMakie.earth(); npoints = 100, z_level = 0, axis = (; type = LScene)); +lp = lines!(a, Point3f.(1:10, 1:10, 110); color = :red, linewidth = 2) +cc = cameracontrols(a.scene) +cc.settings.mouse_translationspeed[] = 0.0 +cc.settings.zoom_shift_lookat[] = false +Makie.update_cam!(a.scene, cc) +p.transformation.transform_func[] = transf +lp.transformation.transform_func[] = transf +f + +# Now, you can scroll around the globe, and the axis will not zoom at all! + + +# Create a white background to obscure some lines you're plotting on the globe: + +random_latlong_points = Point2f.(rand(12) .* 360 .- 180, rand(12) .* 180 .- 90) + +f, a, bg_plot = meshimage(-180..180, -90..90, GeoMakie.earth(); npoints = 100, z_level = -20_000, axis = (; type = LScene)) +bg_plot.transformation.transform_func[] = transf + +# Now, plot some lines on the globe: + +lineplot = lines!(a, GO.segmentize(GO.GeodesicSegments(; max_distance = 100_000), random_latlong_points); color = :black, linewidth = 2) +lineplot.transformation.transform_func[] = transf + +f \ No newline at end of file diff --git a/examples/geodesy_raster_animation.jl b/examples/geodesy_raster_animation.jl new file mode 100644 index 00000000..6b3a24fe --- /dev/null +++ b/examples/geodesy_raster_animation.jl @@ -0,0 +1,52 @@ +# # Geodesy + +# GeoMakie integrates Makie's transformation interface and Geodesy.jl. +# Let's get a Raster and set our data up: +using Makie, GeoMakie, CairoMakie +using Rasters, ArchGDAL, RasterDataSources +using Dates + +# First, load the Raster. +tmax_stack = Rasters.RasterSeries(WorldClim{Climate}, :tmax; month = 1:12) .|> x -> x[:, :, 1] +tmax_raster = cat(tmax_stack...; dims = Ti(1:12)) +# Let's make sure our method works, so we can intitialize plots. +# +# First, we extract a slice of the datacube in time: +current_raster = tmax_raster[Ti(axes(tmax_raster, 3)[1])] +# then, convert it using Makie's `convert_arguments` functionality, which replaces e.g. missing with NaN and straightens axes. +# This is what the plotting recipe will ultimately see. +x, y, z = Makie.convert_arguments(Makie.ContinuousSurface(), current_raster) +# We want any place with missing data to show up as black, so we replace all NaNs with 0. +# For the type of data we're using, and the type of visualization, this is a reasonable assumption. +transform_z = replace(z, NaN => 0.0) + +# Let's now construct the plot. +f_ax_pl, title = with_theme(theme_black()) do + f = Figure() + ax = LScene(f[1, 1]) + pl = meshimage!(ax, current_raster; nan_color=:black) + pl.transformation.transform_func[] = GeoMakie.Geodesy.ECEFfromLLA(GeoMakie.Geodesy.WGS84()) + title = Label(f[begin-1, :], "Title", fontsize = 20, font = :bold, tellwidth = false) + return (Makie.FigureAxisPlot(f, ax, pl), title) +end +f_ax_pl +# Having done all this construction, we set the transformation function: +f_ax_pl.axis.dest = GeoMakie.Geodesy.ECEFfromLLA(GeoMakie.Geodesy.WGS84()) +f +# finally, we can record an animation using the same steps as above: +record(f, "earth_temperature_deformations.mp4", axes(env_transfer_raw_raster, 3); framerate = 30) do i + title.text = Dates.monthname(i) + current_raster = tmax_raster[Ti(i)] + x, y, z = Makie.convert_arguments(Makie.ContinuousSurface(), current_raster) + transform_z = replace(z, NaN => 0.0) # we want a solid Earth in the mesh. + surface_mesh = Makie.surface2mesh(x, y, transform_z .* 1e8) # 1e8 is a scale factor to make deformations look good! + p.input_args[1][] = surface_mesh +end +# ![](earth_temperature_deformations.mp4) +#= +```@cardmeta +Title = "Geodesy" +Description = "The ellipsoidal Earth in 3D" +Cover = fig +``` +=# \ No newline at end of file diff --git a/examples/piracy_at_sea.jl b/examples/piracy_at_sea.jl new file mode 100644 index 00000000..0188b55b --- /dev/null +++ b/examples/piracy_at_sea.jl @@ -0,0 +1,100 @@ +using Makie, GeoMakie + +using GeoJSON, QuackIO, DataFrames +using Downloads + +import GeometryOps as GO, GeoInterface as GI, Proj + +shipping_routes_file = download("https://raw.githubusercontent.com/newzealandpaul/Shipping-Lanes/refs/tags/v1.3.1/data/Shipping_Lanes_v1.geojson", "Shipping_Lanes_v1.geojson") +pirate_attacks_file = download("https://raw.githubusercontent.com/newzealandpaul/Maritime-Pirate-Attacks/refs/heads/main/data/csv/pirate_attacks.csv", "pirate_attacks.csv") + +shipping_routes = GeoJSON.read(shipping_routes_file) +pirate_attacks = QuackIO.read_csv(DataFrame, pirate_attacks_file) + +major = shipping_routes.geometry[findfirst(==("Major"), shipping_routes.Type)] +middle = shipping_routes.geometry[findfirst(==("Middle"), shipping_routes.Type)] +minor = shipping_routes.geometry[findfirst(==("Minor"), shipping_routes.Type)] + +pirate_attacks = DataFrame(CSV.File(pirate_attacks_file; dateformat = "YYYY-MM-DD")) + +using Rasters +using Dates + +pirate_attacks.pre2015 = year.(pirate_attacks.date) .<= 2015 + +groups = groupby(pirate_attacks, :pre2015; sort = true) + +r1 = rasterize(sum, tuple.(groups[1].longitude, groups[1].latitude); fill = 1, res = 0.5, crs = Rasters.EPSG(4326)) +r2 = rasterize(sum, tuple.(groups[2].longitude, groups[2].latitude); fill = 1, res = 0.5, crs = Rasters.EPSG(4326)) + +pts = DimIndices(r1)[(!ismissing).(r1)] +attrs + +dt = DimTable(r1; mergedims = (X, Y)) |> DataFrame |> x -> dropmissing(x, :sum) +dt.XY |> x -> reinterpret(Float64, x) + + + +fig = Figure(; size = (1000, 750)) +ax = GeoAxis(fig[1, 1]; dest = "+proj=wintri") +# ax = Axis(fig[1, 1]; aspect = DataAspect()) +background_plot = meshimage!(ax, -180..180, -90..90, reshape([RGBf(231/255, 255/255, 255/255)], 1, 1)) + +using NaturalEarth +country_plot = poly!(ax, naturalearth("admin_0_countries", 110).geometry; color = RGBf(221/255, 234/255, 214/255), strokewidth = 1, strokecolor = :white) + + +major_plot = lines!(ax, major) +middle_plot = lines!(ax, middle) +minor_plot = lines!(ax, minor) + +major_plot.linewidth = 4 +middle_plot.linewidth = 2 + +minor_plot.visible = false + +major_plot.color = RGBAf(165/255, 223/255, 246/255, 1); +middle_plot.color = RGBAf(165/255, 234/255, 255/255, 1); + + +using Clustering +columns_as_points = permutedims(hcat(groups[1].longitude, groups[1].latitude)) +res = Clustering.kmeans(columns_as_points, 30) +centers = splat(tuple).(eachcol(res.centers)) +markersizes = res.counts + +sp1 = scatter!(ax, centers; markersize = sqrt.(markersizes) .* 1.5, color = color = RGBf(156/255, 96/255, 178/255)) + + +columns_as_points = permutedims(hcat(groups[2].longitude, groups[2].latitude)) +res = Clustering.kmeans(columns_as_points, 110) +centers = splat(tuple).(eachcol(res.centers)) +markersizes = res.counts + +sp2 = scatter!(ax, centers; markersize = sqrt.(markersizes) .* 1.5, color = RGBf(65/255, 183/255, 153/255)) + +sp2.alpha = 0.6 + +# ylims!(ax, -60, 84) + +translate!(sp1, 0, 0, 1) +# ax.aspect[] = 360 #= degrees longitude =# / (84 + 60) #= degrees latitude =# + + +leg_gl = GridLayout(fig[2, 1]) + +leg1 = Legend( + leg_gl[1,1], + [MarkerElement(; marker = sp2.marker, color = sp1.color, markersize) for markersize in sqrt.([5, 15, 30, 60]) .* 1.5], + ["1-10", "11-20", "21-40", "41-80"], + "Pirate Attacks 2015-2020" + ) + +leg2 = Legend( + leg_gl[2,1], + [MarkerElement(; marker = sp2.marker, color = sp2.color, markersize) for markersize in sqrt.([5, 15, 30, 60]) .* 1.5], + ["1-10", "11-20", "21-40", "41-80"], + "Pirate Attacks 2010-2015" +) + +fig \ No newline at end of file diff --git a/examples/specialized/true_icosphere.jl b/examples/specialized/true_icosphere.jl new file mode 100644 index 00000000..02c202ed --- /dev/null +++ b/examples/specialized/true_icosphere.jl @@ -0,0 +1,253 @@ +using GeometryBasics, LinearAlgebra + +using CoordinateTransformations +import GeoInterface as GI + +struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation +end + +function (::UnitCartesianFromGeographic)(geographic_point) + # Longitude is directly translatable to a spherical coordinate + # θ (azimuth) + θ = deg2rad(GI.x(geographic_point)) + # The polar angle is 90 degrees minus the latitude + # ϕ (polar angle) + ϕ = deg2rad(90 - GI.y(geographic_point)) + # Since this is the unit sphere, the radius is assumed to be 1, + # and we don't need to multiply by it. + return Point3( + sin(ϕ) * cos(θ), + sin(ϕ) * sin(θ), + cos(ϕ) + ) +end + +struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation +end + +function (::GeographicFromUnitCartesian)(xyz::AbstractVector) + @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector" + x, y, z = xyz + return Point2( + atan(y, x), + atan(hypot(x, y), z), + ) +end + +struct GeographicUVFromUnitCartesian <: CoordinateTransformations.Transformation +end + +function (::GeographicUVFromUnitCartesian)(xyz::AbstractVector) + @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector" + x, y, z = xyz + return Point2( + ((atan(y, x))/(2pi) + 1)/1.5, + ((atan(hypot(x, y), z))/pi + 0.5)/1.5, + ) +end + + +function add_vertex!(vertex_list, vertex) + push!(vertex_list, LinearAlgebra.normalize(vertex)) + return length(vertex_list) +end + +function get_middle_point!(vertex_list, midpoint_cache, p1::Int, p2::Int) + first_is_smaller = p1 < p2 + smaller_index = first_is_smaller ? p1 : p2 + larger_index = first_is_smaller ? p2 : p1 + key = smaller_index << 32 + larger_index + + if haskey(midpoint_cache, key) + return midpoint_cache[key] + else + midpoint = @. (vertex_list[p1] + vertex_list[p2]) / 2 + midpoint_idx = add_vertex!(vertex_list, midpoint) + midpoint_cache[key] = midpoint_idx + return midpoint_idx + end +end + +function refined_faces!(face, vertex_list, midpoint_cache) + p1, p2, p3 = face + p12 = get_middle_point!(vertex_list, midpoint_cache, p1, p2) + p13 = get_middle_point!(vertex_list, midpoint_cache, p1, p3) + p23 = get_middle_point!(vertex_list, midpoint_cache, p2, p3) + + return ( + TriangleFace{Int}(p1, p12, p13), + TriangleFace{Int}(p2, p23, p12), + TriangleFace{Int}(p3, p13, p23), + TriangleFace{Int}(p12, p23, p13), + ) +end + + +# Helper functions +function crosses_antimeridian(vertices, v1, v2, v3) + s1, s2, s3 = GeographicFromUnitCartesian().(vertices[Vec3(v1, v2, v3)]) + if isapprox(s1[1], pi) || abs.(v1) == Point3{Float64}(0, 0, 1) + return v1 + elseif isapprox(s2[1], pi) || abs.(v2) == Point3{Float64}(0, 0, 1) + return v2 + elseif isapprox(s3[1], pi) || abs.(v3) == Point3{Float64}(0, 0, 1) + return v3 + else + return 0 + end +end + +function split_face_through_vertex(face, antimeridian_vertex, vertices, uvs) + v1, v2, v3 = vertices[face[1]], vertices[face[2]], vertices[face[3]] + # Find the index of the antimeridian vertex in the face + antimeridian_index = findfirst(v -> v == antimeridian_vertex, face) + + if antimeridian_index === nothing + return [face] # No split needed + end + + # Reorder vertices so that the antimeridian vertex is first + ordered_face = circshift(face, 1 - antimeridian_index) + v1, v2, v3 = vertices[ordered_face[1]], vertices[ordered_face[2]], vertices[ordered_face[3]] + uv1, uv2, uv3 = uvs[ordered_face[1]], uvs[ordered_face[2]], uvs[ordered_face[3]] + + # Compute a new vertex along the edge opposite the antimeridian vertex + new_vertex = normalize(@. (v2 + v3)/2) + nv1 = add_vertex!(vertices, new_vertex) + nv2 = add_vertex!(vertices, new_vertex) + + # Calculate UV coordinates for the new vertices + new_uv1 = Vec2{Float64}(0.0, (uv1[2] + uv2[2]) / 2) # Left edge of texture + new_uv2 = Vec2{Float64}(1.0, (uv1[2] + uv3[2]) / 2) # Right edge of texture + push!(uvs, new_uv1) + push!(uvs, new_uv2) + + # Create two new faces, each with a unique set of vertices + new_face1 = (ordered_face[1], ordered_face[2], nv1) # Using new_vertex1 + new_face2 = (ordered_face[1], nv2, ordered_face[3]) # Using new_vertex2 + + return [new_face1, new_face2] +end + +function cartesian_to_spherical(v) + x, y, z = v + θ = acos(z) + ϕ = atan(y, x) + return Point3{Float64}(1, ϕ, θ) +end + +function spherical_to_cartesian(v) + r, ϕ, θ = v + x = r * cos(θ) * cos(ϕ) + y = r * cos(θ) * sin(ϕ) + z = r * cos(θ) + return Point3{Float64}(x, y, z) +end + + + +""" + icosphere(n_refinements::Int = 2) + +Create an icosphere with the specified number of refinements. + +Returns a tuple containing a vector of points (vertices) and a vector of faces. + +## Usage examples +```julia +vertices, faces = icosphere(2) +msh = GeometryBasics.Mesh(vertices, faces) +# you can now plot this however you like: +using GLMakie +wireframe(msh) +``` + +If you want to assign UV coordinates to the vertices, you can compute them from the vertices directly. +In this case we'll map the UVs to a lon-lat equirectangular projection, with `u` going from 0 to 1 +horizontally around the sphere, and `v` going from 0 to 1 vertically from the south pole to the north pole. +```julia +# create the icosphere mesh +vertices, faces = icosphere(2) +# compute the UV coordinates +uvs = [Vec2{Float64}(atan(p[2], p[1]) / (2π) + 0.5, asin(p[3]) / π + 0.5) for p in vertices] + +# create the mesh with UV coordinates +mesh = GeometryBasics.Mesh(GeometryBasics.meta(vertices; uv = uvs), faces) +``` +""" +function icosphere(n_refinements::Int = 2) + # Potential optimizations: + # Precompute the number of vertices and faces based on n_refinements + # Since we know both the input and output number of vertices and faces, + # we can preallocate the arrays using `sizehint!`. + # + # + + ϕ = (1+√5)/2 + vertices = LinearAlgebra.normalize.(Point3{Float64}[ + (-1, ϕ, 0), + ( 1, ϕ, 0), + (-1, -ϕ, 0), + ( 1, -ϕ, 0), + (0, -1, ϕ), + (0, 1, ϕ), + (0, -1, -ϕ), + (0, 1, -ϕ), + ( ϕ, 0, -1), + ( ϕ, 0, 1), + ( -ϕ, 0, -1), + ( -ϕ, 0, 1), + + ]) + + + faces = TriangleFace{Int}.([(1, 12, 6), + (1, 8, 11), + (1, 11, 12), + (2, 6, 10), + (6, 12, 5), + (12, 11, 3), + (11, 8, 7), + (8, 2, 9), + (4, 10, 5), + (4, 7, 9), + (4, 9, 10), + (6, 5, 10), + (3, 5, 12), + (7, 3, 11), + (9, 7, 8), + (10, 9, 2), + ]) + + # Split the faces that aren't already split by the antimeridian + + + midpoint_cache = Dict{Int, Int}() + new_faces = TriangleFace{Int}[] + for _ in 1:n_refinements + for face in faces + append!(new_faces, refined_faces!(face, vertices, midpoint_cache)) + end + tmp = faces + faces = new_faces + new_faces = tmp + empty!(new_faces) + end + + # TODO: add the following corrections to the faces: + # - split all faces along the antimeridian + # - add the ability to refine more points along the poles, to resolve detail. + + # Split faces along the antimeridian (-180 degrees) + uvs = similar(vertices, Vec2f) + antimeridian_faces = TriangleFace{Int}[] + for face in faces + v1, v2, v3 = vertices[face[1]], vertices[face[2]], vertices[face[3]] + + uvs[face[1]] = GeographicUVFromUnitCartesian()(v1) + uvs[face[2]] = GeographicUVFromUnitCartesian()(v2) + uvs[face[3]] = GeographicUVFromUnitCartesian()(v3) + end + + return GeometryBasics.meta(vertices; uv = uvs), faces +end \ No newline at end of file diff --git a/examples/usa_inset.jl b/examples/usa_inset.jl new file mode 100644 index 00000000..ad15541a --- /dev/null +++ b/examples/usa_inset.jl @@ -0,0 +1,94 @@ +#= +# Inset map of the USA + +This example shows how to create an inset map of the USA, that preserves areas as best as possible. + +This example is based on https://geocompx.github.io/geocompkg/articles/us-map.html + +=# + +using GeoMakie, GLMakie +import GeoDataFrames +import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT + +#= +# Data preparation + +The first step is to decide on the best projection for each individual inset. +For this case, we decided to use equal area projections for the maps of the +contiguous 48 states, Hawaii, and Alaska. +While the dataset of Hawaii and Alaska already have this type of projections, +we still need to reproject the `us_states` object to US National Atlas Equal Area: +=# + +using NaturalEarth, DataFrames +us_states = DataFrame(naturalearth("admin_1_states_provinces", 50)) +filter!(:gu_a3 => ==("USA"), us_states) +# Work around https://github.com/JuliaGeo/GeometryOps.jl/issues/215 +us_states = us_states[!, [:geometry, :name]] +# +us_states2163 = GO.reproject(us_states, GFT.EPSG(4326), GFT.EPSG(2163)) +# We also extract the two states of Hawaii and Alaska from the dataset. +hawaii = filter(:name => ==("Hawaii"), us_states2163).geometry |> only +hawaii = GO.reproject(hawaii, GI.crs(hawaii), GFT.ProjString("+proj=aea +lat_0=13 +lon_0=-157 +lat_1=8 +lat_2=18 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs")) + +alaska = filter(:name => ==("Alaska"), us_states2163).geometry |> only +alaska = GO.reproject(alaska, GI.crs(alaska), GFT.EPSG(3467)) + +#= +# Ratio calculation +The second step is to calculate scale relations between the main map +(the contiguous 48 states) and Hawaii, and between the main map and Alaska. +To do so we can calculate areas of the bounding box of each object: +=# + +us_states_range = GI.extent(us_states2163.geometry |> GI.GeometryCollection).Y[2] - GI.extent(us_states2163.geometry |> GI.GeometryCollection).Y[1] +hawaii_range = GI.extent(hawaii).Y[2] - GI.extent(hawaii).Y[1] +alaska_range = GI.extent(alaska).Y[2] - GI.extent(alaska).Y[1] + +# Next, we can calculate the ratio between their areas: + +us_states_hawaii_ratio = hawaii_range / us_states_range +us_states_alaska_ratio = alaska_range / us_states_range +(; us_states_hawaii_ratio, us_states_alaska_ratio) # hide + +#= +# Map creation + +We can now create the inset maps. + +Here, we diverge from the original example, since Makie.jl does not +support creating axes independently of a Figure easily. + +=# + +# First, we instantiate the figure and the axes: +fig = Figure(size = (1200, 800)) +## Alaska takes the top row +ax_alaska = GeoAxis(fig[1, 1]; source = GI.crs(alaska), dest = GI.crs(alaska), tellheight = false, tellwidth = false) +poly!(ax_alaska, alaska; source = GI.crs(alaska), color = :lightgray, strokewidth = 0.75, strokecolor = :darkgray) +## The contiguous 48 states take the bottom row +ax_cont48 = GeoAxis(fig[2, 1]; source = GI.crs(us_states2163), dest = GI.crs(us_states2163), tellheight = false, tellwidth = false) +poly!(ax_cont48, filter(:name => ∉(("Alaska", "Hawaii")), us_states2163).geometry; source = GI.crs(us_states2163), color = :lightgray, strokewidth = 0.75, strokecolor = :darkgray) +## Hawaii will be an inset, so we don't assign it a grid cell yet: +ax_hawaii = GeoAxis(fig; source = GI.crs(hawaii), dest = GI.crs(hawaii), tellheight = false, tellwidth = false) +poly!(ax_hawaii, hawaii; source = GI.crs(hawaii), color = :lightgray, strokewidth = 0.75, strokecolor = :darkgray) + +hidedecorations!(ax_alaska) +hidedecorations!(ax_cont48) +hidedecorations!(ax_hawaii) +fig +# Now, we can set the row heights: +rowsize!(fig.layout, 1, Auto(false, us_states_alaska_ratio)) +rowsize!(fig.layout, 2, Auto(false, 1)) +rowgap!(fig.layout, 0) +fig + +# Now, we move Hawaii to its rightful place: +fig[2, 1] = ax_hawaii +ax_hawaii.valign[] = 0.07 +ax_hawaii.halign[] = 0 +ax_hawaii.height[] = Auto(false, us_states_hawaii_ratio / (us_states_alaska_ratio + 1)) +#= + +=# diff --git a/src/sphere/globeaxis.jl b/src/sphere/globeaxis.jl new file mode 100644 index 00000000..c2de8b50 --- /dev/null +++ b/src/sphere/globeaxis.jl @@ -0,0 +1,356 @@ +#= +# GlobeAxis + +The idea behind the globe axis is to have a 3D axis that can show a globe or a sphere, centered at (0, 0, 0) in cartesian space. + +We can have multiple parametrizations of the globe, from simple unit sphere to fancy ellipsoids or even irregular ellipsoid like things (if you want to plot GRACE data for example). + +Per plot we want the following config: +- source CRS / datum +- default radius / z-value - this can plug into the transformation we give it, so that people can layer their plots + +Per axis we want the following config: +- whether to show default: + - earth map + - smaller mesh sphere (maybe using icosphere) + - coastlines +- whether to show grid lines +- whether to show (potentially) lat/long labels where the grid lines vanish? That would be tricky but a great feature. +- camera controls (rotation speed etc) +- max/min zoom + +For now, a simple implementation is to have a single globe. We can always decide on more or less complexity later. +=# + + +Makie.@Block GlobeAxis <: Makie.AbstractAxis begin + scene::Scene + globe_limits::Observable{Rect3d} + mouseeventhandle::Makie.MouseEventHandle + scrollevents::Observable{Makie.ScrollEvent} + keysevents::Observable{Makie.KeysEvent} + # interactions::Dict{Symbol, Tuple{Bool, Any}} + elements::Dict{Symbol, Any} + final_ellipsoid::Observable{Geodesy.Ellipsoid} + transform_func::Observable{Any} + inv_transform_func::Observable{Any} + @attributes begin + # unused - only for compat with Makie AbstractAxis functions + xscale = identity + yscale = identity + zscale = identity + # Layout observables for Block + "The horizontal alignment of the block in its suggested bounding box." + halign = :center + "The vertical alignment of the block in its suggested bounding box." + valign = :center + "The width setting of the block." + width = Makie.Auto() + "The height setting of the block." + height = Makie.Auto() + "Controls if the parent layout can adjust to this block's width" + tellwidth::Bool = true + "Controls if the parent layout can adjust to this block's height" + tellheight::Bool = true + "The align mode of the block in its parent GridLayout." + alignmode = Makie.Inside() + + # Projection + "Projection of the source data. This is the value plots will default to, but can be overwritten via `plot(...; source=...)`" + source = "+proj=longlat +datum=WGS84" + "Projection that the axis uses to display the data." + dest = "+proj=cart +datum=WGS84" + + + # appearance controls + "The set of fonts which text in the axis should use.s" + fonts = (; regular = "TeX Gyre Heros Makie") + "The axis title string." + title = "" + "The font family of the title." + titlefont = :bold + "The title's font size." + titlesize::Float64 = @inherit(:fontsize, 16f0) + "The gap between axis and title." + titlegap::Float64 = 4f0 + "Controls if the title is visible." + titlevisible::Bool = true + "The horizontal alignment of the title." + titlealign::Symbol = :center + "The color of the title" + titlecolor::RGBAf = @inherit(:textcolor, :black) + "The axis title line height multiplier." + titlelineheight::Float64 = 1 + "The axis subtitle string." + subtitle = "" + "The font family of the subtitle." + subtitlefont = :regular + "The subtitle's font size." + subtitlesize::Float64 = @inherit(:fontsize, 16f0) + "The gap between subtitle and title." + subtitlegap::Float64 = 0 + "Controls if the subtitle is visible." + subtitlevisible::Bool = true + "The color of the subtitle" + subtitlecolor::RGBAf = @inherit(:textcolor, :black) + "The axis subtitle line height multiplier." + subtitlelineheight::Float64 = 1 + + # mesh background + backgroundvisible::Bool = true + "The color of the background mesh." + backgroundcolor = @inherit(:backgroundcolor, :white) + "The opacity of the background mesh." + backgroundalpha::Float64 = 1.0 + "Whether the background mesh has transparency or not." + backgroundtransparency::Bool = true + "The ratio of the background mesh radius to the globe radius." + backgroundradius::Float64 = 0.9 + + # grid lines and ticks + "The x (longitude) ticks - can be a vector or a Makie tick finding algorithm." + xticks = Makie.automatic + "The y (latitude) ticks - can be a vector or a Makie tick finding algorithm." + yticks = Makie.automatic + + "Format for x (longitude) ticks." + xtickformat = Makie.automatic + "Format for y (latitude) ticks." + ytickformat = Makie.automatic + "The font family of the xticklabels." + xticklabelfont = :regular + "The font family of the yticklabels." + yticklabelfont = :regular + "The color of xticklabels." + xticklabelcolor::RGBAf = @inherit(:textcolor, :black) + "The color of yticklabels." + yticklabelcolor::RGBAf = @inherit(:textcolor, :black) + "The font size of the xticklabels." + xticklabelsize::Float64 = @inherit(:fontsize, 16f0) + "The font size of the yticklabels." + yticklabelsize::Float64 = @inherit(:fontsize, 16f0) + "Controls if the xticklabels are visible." + xticklabelsvisible::Bool = true + "Controls if the yticklabels are visible." + yticklabelsvisible::Bool = true + "The space reserved for the xticklabels." + xticklabelspace::Union{Makie.Automatic, Float64} = Makie.automatic + "The space reserved for the yticklabels." + yticklabelspace::Union{Makie.Automatic, Float64} = Makie.automatic + "The space between xticks and xticklabels." + xticklabelpad::Float64 = 5f0 + "The space between yticks and yticklabels." + yticklabelpad::Float64 = 5f0 + "The counterclockwise rotation of the xticklabels in radians." + xticklabelrotation::Float64 = 0f0 + "The counterclockwise rotation of the yticklabels in radians." + yticklabelrotation::Float64 = 0f0 + "The horizontal and vertical alignment of the xticklabels." + xticklabelalign::Union{Makie.Automatic, Tuple{Symbol, Symbol}} = Makie.automatic + "The horizontal and vertical alignment of the yticklabels." + yticklabelalign::Union{Makie.Automatic, Tuple{Symbol, Symbol}} = Makie.automatic + "The size of the xtick marks." + xticksize::Float64 = 6f0 + "The size of the ytick marks." + yticksize::Float64 = 6f0 + "Controls if the xtick marks are visible." + xticksvisible::Bool = true + "Controls if the ytick marks are visible." + yticksvisible::Bool = true + "The alignment of the xtick marks relative to the axis spine (0 = out, 1 = in)." + xtickalign::Float64 = 0f0 + "The alignment of the ytick marks relative to the axis spine (0 = out, 1 = in)." + ytickalign::Float64 = 0f0 + "The width of the xtick marks." + xtickwidth::Float64 = 1f0 + "The width of the ytick marks." + ytickwidth::Float64 = 1f0 + "The color of the xtick marks." + xtickcolor::RGBAf = RGBf(0, 0, 0) + "The color of the ytick marks." + ytickcolor::RGBAf = RGBf(0, 0, 0) + # "The width of the axis spines." + # spinewidth::Float64 = 1f0 + "Controls if the x grid lines are visible." + xgridvisible::Bool = true + "Controls if the y grid lines are visible." + ygridvisible::Bool = true + "The width of the x grid lines." + xgridwidth::Float64 = 1f0 + "The width of the y grid lines." + ygridwidth::Float64 = 1f0 + "The color of the x grid lines." + xgridcolor::RGBAf = RGBAf(0, 0, 0, 0.5) + "The color of the y grid lines." + ygridcolor::RGBAf = RGBAf(0.0, 0, 0, 0.5) + "The linestyle of the x grid lines." + xgridstyle = nothing + "The linestyle of the y grid lines." + ygridstyle = nothing + "Controls if minor ticks on the x axis are visible" + xminorticksvisible::Bool = false + "The alignment of x minor ticks on the axis spine" + xminortickalign::Float64 = 0f0 + "The tick size of x minor ticks" + xminorticksize::Float64 = 4f0 + "The tick width of x minor ticks" + xminortickwidth::Float64 = 1f0 + "The tick color of x minor ticks" + xminortickcolor::RGBAf = :black + "The tick locator for the x minor ticks" + xminorticks = IntervalsBetween(2) + "Controls if minor ticks on the y axis are visible" + yminorticksvisible::Bool = false + "The alignment of y minor ticks on the axis spine" + yminortickalign::Float64 = 0f0 + "The tick size of y minor ticks" + yminorticksize::Float64 = 4f0 + "The tick width of y minor ticks" + yminortickwidth::Float64 = 1f0 + "The tick color of y minor ticks" + yminortickcolor::RGBAf = :black + "The tick locator for the y minor ticks" + yminorticks = IntervalsBetween(2) + "Controls if the x minor grid lines are visible." + xminorgridvisible::Bool = false + "Controls if the y minor grid lines are visible." + yminorgridvisible::Bool = false + "The width of the x minor grid lines." + xminorgridwidth::Float64 = 1f0 + "The width of the y minor grid lines." + yminorgridwidth::Float64 = 1f0 + "The color of the x minor grid lines." + xminorgridcolor::RGBAf = RGBAf(0, 0, 0, 0.05) + "The color of the y minor grid lines." + yminorgridcolor::RGBAf = RGBAf(0, 0, 0, 0.05) + "The linestyle of the x minor grid lines." + xminorgridstyle = nothing + "The linestyle of the y minor grid lines." + yminorgridstyle = nothing + # "Controls if the axis spine is visible." + # spinevisible::Bool = true + # "The color of the axis spine." + # spinecolor::RGBAf = :black + # spinetype::Symbol = :geospine + end +end + +function Makie.initialize_block!(axis::GlobeAxis; scenekw = NamedTuple()) + + # Set up transformations first, so that the scene can be set up + # and linked to those. + transform_obs = Observable{Any}(identity; ignore_equal_values=true) + transform_inv_obs = Observable{Any}(identity; ignore_equal_values=true) + transform_ticks_obs = Observable{Any}(identity; ignore_equal_values=true) + transform_ticks_inv_obs = Observable{Any}(identity; ignore_equal_values=true) + setfield!(axis, :transform_func, transform_obs) + setfield!(axis, :inv_transform_func, transform_inv_obs) + + scene = axis_setup!(axis, transform_obs; scenekw) + + +end + +function reset_limits!(axis::GlobeAxis) + notify(axis.scene.theme.limits) + center!(axis.scene) + return +end +tightlimits!(::GlobeAxis) = nothing # TODO implement!? By getting the bbox of the sphere / ellipsoid and using that to compute the camera eyeposition / lookat / fov + +function axis_setup!(axis::GlobeAxis, transform_obs::Observable; scenekw = NamedTuple()) + + topscene = axis.blockscene + + # pick a camera and draw axis. + scenekw = merge((clear = false, camera=cam3d!), scenekw) + axis.scene = Scene(blockscene, lift(round_to_IRect2D, blockscene, axis.layoutobservables.computedbbox); transformation = Makie.Transformation(topscene.transformation; transform_func = transform_obs), scenekw...) + # Axis should not have a transform func, but all other plots should! + _ax = axis.scene[OldAxis] + isnothing(_ax) || (_ax.transformation = Makie.Transformation(axis.scene; transform_func = identity)) + + on(blockscene, axis.show_axis) do show_axis + ax = axis.scene[OldAxis] + if show_axis + if isnothing(ax) + # Add axis on first plot!, if requested + # update limits when scene limits change + limits = lift(blockscene, axis.scene.theme.limits) do lims + if lims === automatic + dl = boundingbox(axis.scene, p -> Makie.isaxis(p) || Makie.not_in_data_space(p)) + if any(isinf, widths(dl)) || any(isinf, Makie.origin(dl)) + Rect3d((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)) + else + dl + end + else + lims + end + end + Makie.axis3d!(axis.scene, limits) + # Make sure axis is always in pos 1 + sort!(axis.scene.plots, by=!Makie.isaxis) + else + ax.visible = true + end + else + if !isnothing(ax) + ax.visible = false + end + end + end + notify(axis.show_axis) + + return axis.scene +end + + +# This is where we override the stuff to make it our stuff. +function Makie.plot!(axis::GeoAxis, plot::Makie.AbstractPlot) + source = pop!(plot.kw, :source, axis.source) + transformfunc = lift(create_transform, axis.dest, source) + trans = Transformation(transformfunc; get(plot.kw, :transformation, Attributes())...) + plot.kw[:transformation] = trans + Makie.plot!(axis.scene, plot) + # some area-like plots basically always look better if they cover the whole plot area. + # adjust the limit margins in those cases automatically. + # However, for spheres, we want to keep user zoom level if possible. + # Makie.needs_tight_limits(plot) && Makie.tightlimits!(axis) + if Makie.is_open_or_any_parent(axis.scene) + Makie.reset_limits!(axis) + end + return plot +end +function Makie.MakieCore._create_plot!(F, attributes::Dict, ax::GeoAxis, args...) + source = pop!(attributes, :source, nothing) + dest = pop!(attributes, :dest, nothing) + plot = Plot{Makie.default_plot_func(F, args)}(args, attributes) + isnothing(source) || (plot.kw[:source] = source) + isnothing(dest) || (plot.kw[:dest] = dest) + Makie.plot!(ax, plot) + return plot +end + +# Makie.interactions(ax::GlobeAxis) = ax.interactions + +function Makie.update_state_before_display!(ax::GlobeAxis) + Makie.reset_limits!(ax) + return +end + + +# Legend API + +function Makie.get_plots(ax::GlobeAxis) + n_skip = isnothing(ax.scene[OldAxis]) ? 1 : 2 + return Makie.get_plots(ax.scene)[n_skip:end] +end + +function Makie.Legend(fig_or_scene, axis::GlobeAxis, title = nothing; merge = false, unique = false, kwargs...) + plots, labels = Makie.get_labeled_plots(axis, merge = merge, unique = unique) + isempty(plots) && error("There are no plots with labels in the given axis that can be put in the legend. Supply labels to plotting functions like `plot(args...; label = \"My label\")`") + Makie.Legend(fig_or_scene, plots, labels, title; kwargs...) +end + + + diff --git a/src/sphere/icosphere.jl b/src/sphere/icosphere.jl new file mode 100644 index 00000000..ea0c3fd3 --- /dev/null +++ b/src/sphere/icosphere.jl @@ -0,0 +1,210 @@ +using GeometryBasics, LinearAlgebra + +using CoordinateTransformations +import GeoInterface as GI + + +function add_vertex!(vertex_list, vertex) + push!(vertex_list, LinearAlgebra.normalize(vertex)) + return length(vertex_list) +end + +function get_middle_point!(vertex_list, midpoint_cache, p1::Int, p2::Int) + first_is_smaller = p1 < p2 + smaller_index = first_is_smaller ? p1 : p2 + larger_index = first_is_smaller ? p2 : p1 + key = smaller_index << 32 + larger_index + + if haskey(midpoint_cache, key) + return midpoint_cache[key] + else + midpoint = @. (vertex_list[p1] + vertex_list[p2]) / 2 + midpoint_idx = add_vertex!(vertex_list, midpoint) + midpoint_cache[key] = midpoint_idx + return midpoint_idx + end +end + +function refined_faces!(face, vertex_list, midpoint_cache) + p1, p2, p3 = face + p12 = get_middle_point!(vertex_list, midpoint_cache, p1, p2) + p13 = get_middle_point!(vertex_list, midpoint_cache, p1, p3) + p23 = get_middle_point!(vertex_list, midpoint_cache, p2, p3) + + return ( + TriangleFace{Int}(p1, p12, p13), + TriangleFace{Int}(p2, p23, p12), + TriangleFace{Int}(p3, p13, p23), + TriangleFace{Int}(p12, p23, p13), + ) +end + + +# Helper functions +function crosses_antimeridian(vertices, v1, v2, v3) + s1, s2, s3 = GeographicFromUnitCartesian().(vertices[Vec3(v1, v2, v3)]) + if isapprox(s1[1], pi) || abs.(v1) == Point3{Float64}(0, 0, 1) + return v1 + elseif isapprox(s2[1], pi) || abs.(v2) == Point3{Float64}(0, 0, 1) + return v2 + elseif isapprox(s3[1], pi) || abs.(v3) == Point3{Float64}(0, 0, 1) + return v3 + else + return 0 + end +end + +function split_face_through_vertex(face, antimeridian_vertex, vertices, uvs) + v1, v2, v3 = vertices[face[1]], vertices[face[2]], vertices[face[3]] + # Find the index of the antimeridian vertex in the face + antimeridian_index = findfirst(v -> v == antimeridian_vertex, face) + + if antimeridian_index === nothing + return [face] # No split needed + end + + # Reorder vertices so that the antimeridian vertex is first + ordered_face = circshift(face, 1 - antimeridian_index) + v1, v2, v3 = vertices[ordered_face[1]], vertices[ordered_face[2]], vertices[ordered_face[3]] + uv1, uv2, uv3 = uvs[ordered_face[1]], uvs[ordered_face[2]], uvs[ordered_face[3]] + + # Compute a new vertex along the edge opposite the antimeridian vertex + new_vertex = normalize(@. (v2 + v3)/2) + nv1 = add_vertex!(vertices, new_vertex) + nv2 = add_vertex!(vertices, new_vertex) + + # Calculate UV coordinates for the new vertices + new_uv1 = Vec2{Float64}(0.0, (uv1[2] + uv2[2]) / 2) # Left edge of texture + new_uv2 = Vec2{Float64}(1.0, (uv1[2] + uv3[2]) / 2) # Right edge of texture + push!(uvs, new_uv1) + push!(uvs, new_uv2) + + # Create two new faces, each with a unique set of vertices + new_face1 = (ordered_face[1], ordered_face[2], nv1) # Using new_vertex1 + new_face2 = (ordered_face[1], nv2, ordered_face[3]) # Using new_vertex2 + + return [new_face1, new_face2] +end + +function cartesian_to_spherical(v) + x, y, z = v + θ = acos(z) + ϕ = atan(y, x) + return Point3{Float64}(1, ϕ, θ) +end + +function spherical_to_cartesian(v) + r, ϕ, θ = v + x = r * cos(θ) * cos(ϕ) + y = r * cos(θ) * sin(ϕ) + z = r * cos(θ) + return Point3{Float64}(x, y, z) +end + + + +""" + icosphere(n_refinements::Int = 2) + +Create an icosphere with the specified number of refinements. + +Returns a tuple containing a vector of points (vertices) and a vector of faces. + +## Usage examples +```julia +vertices, faces = icosphere(2) +msh = GeometryBasics.Mesh(vertices, faces) +# you can now plot this however you like: +using GLMakie +wireframe(msh) +``` + +If you want to assign UV coordinates to the vertices, you can compute them from the vertices directly. +In this case we'll map the UVs to a lon-lat equirectangular projection, with `u` going from 0 to 1 +horizontally around the sphere, and `v` going from 0 to 1 vertically from the south pole to the north pole. +```julia +# create the icosphere mesh +vertices, faces = icosphere(2) +# compute the UV coordinates +uvs = [Vec2{Float64}(atan(p[2], p[1]) / (2π) + 0.5, asin(p[3]) / π + 0.5) for p in vertices] + +# create the mesh with UV coordinates +mesh = GeometryBasics.Mesh(GeometryBasics.meta(vertices; uv = uvs), faces) +``` +""" +function icosphere(n_refinements::Int = 2) + # Potential optimizations: + # Precompute the number of vertices and faces based on n_refinements + # Since we know both the input and output number of vertices and faces, + # we can preallocate the arrays using `sizehint!`. + # + # + + ϕ = (1+√5)/2 + vertices = LinearAlgebra.normalize.(Point3{Float64}[ + (-1, ϕ, 0), + ( 1, ϕ, 0), + (-1, -ϕ, 0), + ( 1, -ϕ, 0), + (0, -1, ϕ), + (0, 1, ϕ), + (0, -1, -ϕ), + (0, 1, -ϕ), + ( ϕ, 0, -1), + ( ϕ, 0, 1), + ( -ϕ, 0, -1), + ( -ϕ, 0, 1), + + ]) + + + faces = TriangleFace{Int}.([(1, 12, 6), + (1, 8, 11), + (1, 11, 12), + (2, 6, 10), + (6, 12, 5), + (12, 11, 3), + (11, 8, 7), + (8, 2, 9), + (4, 10, 5), + (4, 7, 9), + (4, 9, 10), + (6, 5, 10), + (3, 5, 12), + (7, 3, 11), + (9, 7, 8), + (10, 9, 2), + ]) + + # Split the faces that aren't already split by the antimeridian + + + midpoint_cache = Dict{Int, Int}() + new_faces = TriangleFace{Int}[] + for _ in 1:n_refinements + for face in faces + append!(new_faces, refined_faces!(face, vertices, midpoint_cache)) + end + tmp = faces + faces = new_faces + new_faces = tmp + empty!(new_faces) + end + + # TODO: add the following corrections to the faces: + # - split all faces along the antimeridian + # - add the ability to refine more points along the poles, to resolve detail. + + # Split faces along the antimeridian (-180 degrees) + uvs = similar(vertices, Vec2f) + antimeridian_faces = TriangleFace{Int}[] + for face in faces + v1, v2, v3 = vertices[face[1]], vertices[face[2]], vertices[face[3]] + + uvs[face[1]] = GeographicUVFromUnitCartesian()(v1) + uvs[face[2]] = GeographicUVFromUnitCartesian()(v2) + uvs[face[3]] = GeographicUVFromUnitCartesian()(v3) + end + + return GeometryBasics.meta(vertices; uv = uvs), faces +end \ No newline at end of file diff --git a/src/sphere/unit_sphere_transforms.jl b/src/sphere/unit_sphere_transforms.jl new file mode 100644 index 00000000..688434ba --- /dev/null +++ b/src/sphere/unit_sphere_transforms.jl @@ -0,0 +1,55 @@ +import CoordinateTransformations +import GeoInterface as GI + +struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation +end + +function (::UnitCartesianFromGeographic)(geographic_point) + # Longitude is directly translatable to a spherical coordinate + # θ (azimuth) + θ = deg2rad(GI.x(geographic_point)) + # The polar angle is 90 degrees minus the latitude, so we must use cosine instead of sine. + # ϕ (polar angle) + ϕ = deg2rad(90 - GI.y(geographic_point)) + # Since this is the unit sphere, the radius is assumed to be 1, + # and we don't need to multiply by it. + return Point3( + sin(ϕ) * cos(θ), + sin(ϕ) * sin(θ), + cos(ϕ) + ) +end + +Base.inv(::UnitCartesianFromGeographic) = GeographicFromUnitCartesian() + + + + + +struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation +end + +function (::GeographicFromUnitCartesian)(xyz::AbstractVector) + @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector" + x, y, z = xyz + return Point2( + atan(y, x), + atan(hypot(x, y), z), + ) +end + +Base.inv(::GeographicFromUnitCartesian) = UnitCartesianFromGeographic() + + + +struct GeographicUVFromUnitCartesian <: CoordinateTransformations.Transformation +end + +function (::GeographicUVFromUnitCartesian)(xyz::AbstractVector) + @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector" + x, y, z = xyz + return Point2( + ((atan(y, x))/(2pi) + 1)/1.5, + ((atan(hypot(x, y), z))/pi + 0.5)/1.5, + ) +end From 2be2e3e43a3ab5c895d0fb8f154a101cc5e54906 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 29 Oct 2024 17:40:58 -0700 Subject: [PATCH 02/12] Improvements to geodesy, add CoordinateTransformations to Project --- Project.toml | 2 ++ examples/geodesy.jl | 52 --------------------------------- examples/specialized/earth3d.jl | 13 --------- src/geodesy.jl | 5 +--- 4 files changed, 3 insertions(+), 69 deletions(-) delete mode 100644 examples/geodesy.jl delete mode 100644 examples/specialized/earth3d.jl diff --git a/Project.toml b/Project.toml index 77081df6..c1f27545 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.7.9" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" @@ -24,6 +25,7 @@ StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" [compat] Colors = "0.12" +CoordinateTransformations = "0.6.3" Downloads = "1" GeoFormatTypes = "0.4" GeoInterface = "0.5, 1.0" diff --git a/examples/geodesy.jl b/examples/geodesy.jl deleted file mode 100644 index 6b3a24fe..00000000 --- a/examples/geodesy.jl +++ /dev/null @@ -1,52 +0,0 @@ -# # Geodesy - -# GeoMakie integrates Makie's transformation interface and Geodesy.jl. -# Let's get a Raster and set our data up: -using Makie, GeoMakie, CairoMakie -using Rasters, ArchGDAL, RasterDataSources -using Dates - -# First, load the Raster. -tmax_stack = Rasters.RasterSeries(WorldClim{Climate}, :tmax; month = 1:12) .|> x -> x[:, :, 1] -tmax_raster = cat(tmax_stack...; dims = Ti(1:12)) -# Let's make sure our method works, so we can intitialize plots. -# -# First, we extract a slice of the datacube in time: -current_raster = tmax_raster[Ti(axes(tmax_raster, 3)[1])] -# then, convert it using Makie's `convert_arguments` functionality, which replaces e.g. missing with NaN and straightens axes. -# This is what the plotting recipe will ultimately see. -x, y, z = Makie.convert_arguments(Makie.ContinuousSurface(), current_raster) -# We want any place with missing data to show up as black, so we replace all NaNs with 0. -# For the type of data we're using, and the type of visualization, this is a reasonable assumption. -transform_z = replace(z, NaN => 0.0) - -# Let's now construct the plot. -f_ax_pl, title = with_theme(theme_black()) do - f = Figure() - ax = LScene(f[1, 1]) - pl = meshimage!(ax, current_raster; nan_color=:black) - pl.transformation.transform_func[] = GeoMakie.Geodesy.ECEFfromLLA(GeoMakie.Geodesy.WGS84()) - title = Label(f[begin-1, :], "Title", fontsize = 20, font = :bold, tellwidth = false) - return (Makie.FigureAxisPlot(f, ax, pl), title) -end -f_ax_pl -# Having done all this construction, we set the transformation function: -f_ax_pl.axis.dest = GeoMakie.Geodesy.ECEFfromLLA(GeoMakie.Geodesy.WGS84()) -f -# finally, we can record an animation using the same steps as above: -record(f, "earth_temperature_deformations.mp4", axes(env_transfer_raw_raster, 3); framerate = 30) do i - title.text = Dates.monthname(i) - current_raster = tmax_raster[Ti(i)] - x, y, z = Makie.convert_arguments(Makie.ContinuousSurface(), current_raster) - transform_z = replace(z, NaN => 0.0) # we want a solid Earth in the mesh. - surface_mesh = Makie.surface2mesh(x, y, transform_z .* 1e8) # 1e8 is a scale factor to make deformations look good! - p.input_args[1][] = surface_mesh -end -# ![](earth_temperature_deformations.mp4) -#= -```@cardmeta -Title = "Geodesy" -Description = "The ellipsoidal Earth in 3D" -Cover = fig -``` -=# \ No newline at end of file diff --git a/examples/specialized/earth3d.jl b/examples/specialized/earth3d.jl deleted file mode 100644 index 5fad41b4..00000000 --- a/examples/specialized/earth3d.jl +++ /dev/null @@ -1,13 +0,0 @@ -using GLMakie, GeoMakie, Geodesy - -transf = Geodesy.ECEFfromLLA(Geodesy.WGS84()) - -f, a, p = meshimage(-180..180, -90..90, GeoMakie.earth(); npoints = 100, z_level = 0, axis = (; type = LScene)); -lp = lines!(a, Point3f.(1:10, 1:10, 110); color = :red, linewidth = 2) -cc = cameracontrols(a.scene) -cc.settings.mouse_translationspeed[] = 0.0 -cc.settings.zoom_shift_lookat[] = false -Makie.update_cam!(a.scene, cc) -p.transformation.transform_func[] = transf -lp.transformation.transform_func[] = transf -f diff --git a/src/geodesy.jl b/src/geodesy.jl index ba0c0245..5adda763 100644 --- a/src/geodesy.jl +++ b/src/geodesy.jl @@ -13,9 +13,6 @@ function Makie.apply_transform(f::Geodesy.ECEFfromLLA, pt::V) where V <: VecType # However, the `LLA` coordinate space expects x to be latitude and y to be longitude, # so we have to manually swap the coordinates. return V((f(LLA(pt[2], pt[1], pt[3])))...) - # Enterprising observers will note the division by 2,500. This is to "normalize" the sphere; - # Geodesy.jl outputs in meters, which would make the sphere hilariously large. - # This also fits well with Makie's limit finding, which works in input space, and not transformed space. end # If a Point2f is passed, we decide to handle that by assuming altitude to be 0. @@ -36,7 +33,7 @@ end Makie.inverse_transform(f::Geodesy.ECEFfromLLA) = Base.inv(f) # and its application: function Makie.apply_transform(f::Geodesy.LLAfromECEF, pt::V) where V <: VecTypes{3, T} where {T} - return V((f(ECEF(pt[1], pt[2], pt[3])))...) # invert the previous scale factor + return V((f(ECEF(pt[1], pt[2], pt[3])))...) end function Makie.apply_transform(f::Geodesy.LLAfromECEF, pt::V) where V <: VecTypes{N, T} where {N, T} From b0bc40639a08070cfc1c2b622bb85a168cfbd6d1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 20 Nov 2024 16:31:22 -0500 Subject: [PATCH 03/12] Create "globe transforms" that guarantee 2d->3d TODO: set poly_convert to subsample here. --- src/sphere/globetransform.jl | 145 +++++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) create mode 100644 src/sphere/globetransform.jl diff --git a/src/sphere/globetransform.jl b/src/sphere/globetransform.jl new file mode 100644 index 00000000..5882a4a7 --- /dev/null +++ b/src/sphere/globetransform.jl @@ -0,0 +1,145 @@ +""" + abstract type GlobeTransform + +The supertype for all globe transforms. All subtypes must be callable with a Point, and have a field `zlevel`. + +Construct via [`create_globe_transform`](@ref). Used in [`GlobeAxis`](@ref). +""" +abstract type GlobeTransform end + +function create_globe_transform(dest::Geodesy.Ellipsoid, src::Type, zlevel) + return GeodesyGlobeTransform(dest, src, zlevel) +end + +function create_globe_transform(dest::Geodesy.Ellipsoid, src::Union{String, GeoFormatTypes.GeoFormat}, zlevel) + return ProjGlobeTransform(dest, src, zlevel) +end + +function create_globe_transform(dest, src, zlevel) + return ProjGlobeTransform(dest, src, zlevel) +end + + +struct ProjGlobeTransform <: GlobeTransform + transf::Proj.Transformation + zlevel::Float64 +end + +function ProjGlobeTransform(dest::Geodesy.Datum, src, zlevel) + return ProjGlobeTransform(Geodesy.ellipsoid(dest), src, zlevel) +end + +function ProjGlobeTransform(dest::Geodesy.Ellipsoid, src, zlevel) + dest_str = "+proj=cart +a=$(dest.a) +f=$(dest.f)" + transf = create_transform(dest_str, src) + return ProjGlobeTransform(transf, zlevel) +end + +ProjGlobeTransform(dest, src, zlevel) = ProjGlobeTransform(create_transform(dest, src), zlevel) + +function (f::ProjGlobeTransform)(pt::Makie.VecTypes{3, T}) where {T} + return Point3(f.transf((pt[1], pt[2], pt[3] + f.zlevel))) +end + +function (f::ProjGlobeTransform)(pt::Makie.VecTypes{2, T}) where {T} + return Point3(f.transf((pt[1], pt[2], f.zlevel))) +end + +function (f::ProjGlobeTransform)(pt) + z = GI.is3d(pt) ? GI.z(pt) + f.zlevel : f.zlevel + return Point3(f.transf((GI.x(pt), GI.y(pt), z))) +end + + +struct GeodesyGlobeTransform{CoordType, T} <: GlobeTransform + transf::T + zlevel::Float64 +end + + +function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.LLA}) + return GeodesyGlobeTransform{src}(Geodesy.ECEFfromLLA(dest), zlevel) +end + + +function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.UTM}) + return GeodesyGlobeTransform{src}(Geodesy.ECEFfromUTM(dest), zlevel) +end + +function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.UTMZ}) + return GeodesyGlobeTransform{src}(Geodesy.ECEFfromUTMZ(dest), zlevel) +end + +function (f::GeodesyGlobeTransform{Geodesy.LLA})(pt::Makie.VecTypes{3, T}) where {T} + return Point3(f.transf(Geodesy.LLA(pt[2], pt[1], pt[3] + f.zlevel))...) +end + + +function (f::GeodesyGlobeTransform{Geodesy.LLA})(pt::Makie.VecTypes{2, T}) where {T} + return Point3(f.transf(Geodesy.LLA(pt[2], pt[1], f.zlevel))...) +end + + + +# Generic methods for globe transforms + +function Makie.apply_transform(t::GlobeTransform, v::V) where V <: VecTypes{2,T} where {T} + return t(Point2(v[1], v[2])) +end + + +function Makie.apply_transform(t::GlobeTransform, v::V) where V <: VecTypes{3,T} where {T} + return t(Point3(v[1], v[2], v[3])) +end + + +function iterated_bounds_3d(f, (xmin, xmax), (ymin, ymax), N = 21; zlevel = 0) + umin, umax = Inf, -Inf + vmin, vmax = Inf, -Inf + wmin, wmax = Inf, -Inf + for x in LinRange(xmin, xmax, N) + for y in LinRange(ymin, ymax, N) + u, v, w = Makie.apply_transform(f, Vec(x, y, zlevel)) + isfinite(u) && (umin = min(umin, u)) + isfinite(u) && (umax = max(umax, u)) + isfinite(v) && (vmin = min(vmin, v)) + isfinite(v) && (vmax = max(vmax, v)) + isfinite(w) && (wmin = min(wmin, w)) + isfinite(w) && (wmax = max(wmax, w)) + end + end + return (umin, umax), (vmin, vmax), (wmin, wmax) +end + +# Converting rectangles requires densifying the edges. +function Makie.apply_transform(f::GlobeTransform, r::Rect2{T}; zlevel = f.zlevel) where {T} + xmin, ymin = minimum(r) + xmax, ymax = maximum(r) + + if isapprox(xmin, -180, rtol = 1e-4) + xmin - -180e0 + end + if isapprox(xmax, 180; rtol = 1e-4) + xmax = 180e0 + end + if isapprox(ymin, -90; rtol = 1e-4) + ymin = -90e0 + end + if isapprox(ymax, 90; rtol = 1e-4) + ymax = 90e0 + end + + try + (umin, umax), (vmin, vmax), (wmin, wmax) = iterated_bounds_3d(f, (xmin,xmax), (ymin,ymax); zlevel) + return Rect3(Vec3(T(umin), T(vmin), T(wmin)), Vec3(T(umax-umin), T(vmax-vmin), T(wmax-wmin))) + catch e + @show r + rethrow(e) + end +end + +function Makie.apply_transform(f::GlobeTransform, r::Rect3{T}) where {T} + r2 = Rect2{T}(r) + tr2 = Makie.apply_transform(f, r2; zlevel = origin(r)[3] + widths(r)[3] + f.zlevel) + return Rect3{T}(tr2.origin, tr2.widths) +end From 2627c5aa654ffaa3c8e1c65a6ba417d881f6e1d8 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 20 Nov 2024 16:31:35 -0500 Subject: [PATCH 04/12] get GlobeAxis working --- src/GeoMakie.jl | 5 ++ src/sphere/globeaxis.jl | 158 ++++++++++++++++++++++++---------------- test/globeaxis.jl | 2 + test/globetransforms.jl | 0 4 files changed, 103 insertions(+), 62 deletions(-) create mode 100644 test/globeaxis.jl create mode 100644 test/globetransforms.jl diff --git a/src/GeoMakie.jl b/src/GeoMakie.jl index 8a520cdc..77b31bbe 100644 --- a/src/GeoMakie.jl +++ b/src/GeoMakie.jl @@ -60,6 +60,11 @@ include("makie-axis.jl") include("mesh_image.jl") include("linesplitting.jl") +include("sphere/unit_sphere_transforms.jl") +include("sphere/icosphere.jl") +include("sphere/globetransform.jl") +include("sphere/globeaxis.jl") + include("triangulation3d.jl") @reexport using Colors, Makie diff --git a/src/sphere/globeaxis.jl b/src/sphere/globeaxis.jl index c2de8b50..139b88a9 100644 --- a/src/sphere/globeaxis.jl +++ b/src/sphere/globeaxis.jl @@ -24,14 +24,16 @@ For now, a simple implementation is to have a single globe. We can always decid Makie.@Block GlobeAxis <: Makie.AbstractAxis begin + @forwarded_layout scene::Scene + lscene::LScene globe_limits::Observable{Rect3d} - mouseeventhandle::Makie.MouseEventHandle - scrollevents::Observable{Makie.ScrollEvent} - keysevents::Observable{Makie.KeysEvent} + # mouseeventhandle::Makie.MouseEventHandle + # scrollevents::Observable{Makie.ScrollEvent} + # keysevents::Observable{Makie.KeysEvent} # interactions::Dict{Symbol, Tuple{Bool, Any}} elements::Dict{Symbol, Any} - final_ellipsoid::Observable{Geodesy.Ellipsoid} + ellipsoid::Observable{Geodesy.Ellipsoid} transform_func::Observable{Any} inv_transform_func::Observable{Any} @attributes begin @@ -63,6 +65,8 @@ Makie.@Block GlobeAxis <: Makie.AbstractAxis begin # appearance controls + "Controls if the axis is visible. This is a regular 3D OldAxis for now." + show_axis::Bool = true "The set of fonts which text in the axis should use.s" fonts = (; regular = "TeX Gyre Heros Makie") "The axis title string." @@ -235,8 +239,15 @@ Makie.@Block GlobeAxis <: Makie.AbstractAxis begin end end -function Makie.initialize_block!(axis::GlobeAxis; scenekw = NamedTuple()) +_geodesy_ellipsoid_from(g::Geodesy.Ellipsoid) = g +_geodesy_ellipsoid_from(g::Geodesy.Datum) = Geodesy.ellipsoid(g) +_geodesy_ellipsoid_from(n::NamedTuple{(:a, :f), Tuple{Float64, Float64}}) = Geodesy.ellipsoid(; n...) +function Makie.initialize_block!(axis::GlobeAxis; ellipsoid = Geodesy.wgs84_ellipsoid, scenekw = NamedTuple()) + + # axis.show_axis = true + # TODO: set up globe ellipsoid ahead of time. + setfield!(axis, :ellipsoid, Observable{Geodesy.Ellipsoid}(_geodesy_ellipsoid_from(ellipsoid))) # TODO: change this to Geodesy.Datum, allow for a Geodesy.FlexibleDatum type that takes an ellipsoid. # Set up transformations first, so that the scene can be set up # and linked to those. transform_obs = Observable{Any}(identity; ignore_equal_values=true) @@ -246,72 +257,93 @@ function Makie.initialize_block!(axis::GlobeAxis; scenekw = NamedTuple()) setfield!(axis, :transform_func, transform_obs) setfield!(axis, :inv_transform_func, transform_inv_obs) - scene = axis_setup!(axis, transform_obs; scenekw) + lscene = LScene(axis.layout[1, 1]; show_axis = axis.show_axis, scenekw) + scene = lscene.scene + setfield!(axis, :scene, scene) + setfield!(axis, :lscene, lscene) + + # set up the camera for the axis + cc = cameracontrols(scene) + cc.settings.mouse_translationspeed[] = 0.0 + cc.settings.zoom_shift_lookat[] = false + cc.lookat[] = Makie.Vec3f(0, 0, 0) # center the camera on the globe + # TODO: decide - should we use cam3d_cad or oldcam3d here? + Makie.update_cam!(scene, cc) + setfield!(axis, :globe_limits, Observable{Rect3d}(Makie.boundingbox(axis.scene))) + setfield!(axis, :elements, Dict{Symbol, Any}()) + + return end -function reset_limits!(axis::GlobeAxis) - notify(axis.scene.theme.limits) - center!(axis.scene) +function Makie.reset_limits!(axis::GlobeAxis; kwargs...) + Makie.reset_limits!(axis.lscene; kwargs...) return end tightlimits!(::GlobeAxis) = nothing # TODO implement!? By getting the bbox of the sphere / ellipsoid and using that to compute the camera eyeposition / lookat / fov -function axis_setup!(axis::GlobeAxis, transform_obs::Observable; scenekw = NamedTuple()) - - topscene = axis.blockscene - - # pick a camera and draw axis. - scenekw = merge((clear = false, camera=cam3d!), scenekw) - axis.scene = Scene(blockscene, lift(round_to_IRect2D, blockscene, axis.layoutobservables.computedbbox); transformation = Makie.Transformation(topscene.transformation; transform_func = transform_obs), scenekw...) - # Axis should not have a transform func, but all other plots should! - _ax = axis.scene[OldAxis] - isnothing(_ax) || (_ax.transformation = Makie.Transformation(axis.scene; transform_func = identity)) - - on(blockscene, axis.show_axis) do show_axis - ax = axis.scene[OldAxis] - if show_axis - if isnothing(ax) - # Add axis on first plot!, if requested - # update limits when scene limits change - limits = lift(blockscene, axis.scene.theme.limits) do lims - if lims === automatic - dl = boundingbox(axis.scene, p -> Makie.isaxis(p) || Makie.not_in_data_space(p)) - if any(isinf, widths(dl)) || any(isinf, Makie.origin(dl)) - Rect3d((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)) - else - dl - end - else - lims - end - end - Makie.axis3d!(axis.scene, limits) - # Make sure axis is always in pos 1 - sort!(axis.scene.plots, by=!Makie.isaxis) - else - ax.visible = true - end - else - if !isnothing(ax) - ax.visible = false - end - end - end - notify(axis.show_axis) - - return axis.scene -end +# function axis_setup!(axis::GlobeAxis, transform_obs::Observable; scenekw = NamedTuple()) + +# topscene = axis.blockscene + +# # pick a camera and draw axis. +# scenekw = merge((clear = false, camera=Makie.cam3d!), scenekw) +# axis.scene = Scene(topscene, lift(Makie.round_to_IRect2D, topscene, axis.layoutobservables.computedbbox); transformation = Makie.Transformation(topscene.transformation; transform_func = transform_obs), scenekw...) +# # Axis should not have a transform func, but all other plots should! +# _ax = axis.scene[Makie.OldAxis] +# isnothing(_ax) || (_ax.transformation = Makie.Transformation(axis.scene; transform_func = identity)) + +# on(topscene, axis.show_axis) do show_axis +# ax = axis.scene[Makie.OldAxis] +# if show_axis +# if isnothing(ax) +# # Add axis on first plot!, if requested +# # update limits when scene limits change +# limits = lift(topscene, axis.scene.theme.limits) do lims +# if lims === Makie.automatic +# dl = boundingbox(axis.scene, p -> Makie.isaxis(p) || Makie.not_in_data_space(p)) +# if any(isinf, widths(dl)) || any(isinf, Makie.origin(dl)) +# Rect3d((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)) +# else +# dl +# end +# else +# lims +# end +# end +# Makie.axis3d!(axis.scene, limits) +# # Make sure axis is always in pos 1 +# sort!(axis.scene.plots, by=!Makie.isaxis) +# else +# ax.visible = true +# end +# else +# if !isnothing(ax) +# ax.visible = false +# end +# end +# end +# notify(axis.show_axis) + +# return axis.scene +# end # This is where we override the stuff to make it our stuff. -function Makie.plot!(axis::GeoAxis, plot::Makie.AbstractPlot) +function Makie.plot!(axis::GlobeAxis, plot::Makie.AbstractPlot) source = pop!(plot.kw, :source, axis.source) - transformfunc = lift(create_transform, axis.dest, source) - trans = Transformation(transformfunc; get(plot.kw, :transformation, Attributes())...) + zlevel = pop!(plot.kw, :zlevel, 0) + # @show plot.kw + transformfunc = lift(create_globe_transform, axis.ellipsoid, source, zlevel) + trans = Makie.Transformation(transformfunc; get(plot.kw, :transformation, Attributes())...) plot.kw[:transformation] = trans + Makie.plot!(axis.scene, plot) + + # reassign popped observables back to attributes, so that they can be called by the user... + # plot.attributes.attributes[:source] = source + # plot.attributes.attributes[:zlevel] = zlevel # some area-like plots basically always look better if they cover the whole plot area. # adjust the limit margins in those cases automatically. # However, for spheres, we want to keep user zoom level if possible. @@ -321,12 +353,14 @@ function Makie.plot!(axis::GeoAxis, plot::Makie.AbstractPlot) end return plot end -function Makie.MakieCore._create_plot!(F, attributes::Dict, ax::GeoAxis, args...) - source = pop!(attributes, :source, nothing) - dest = pop!(attributes, :dest, nothing) +function Makie.MakieCore._create_plot!(F, attributes::Dict, ax::GlobeAxis, args...) + source = pop!(attributes, :source, ax.source) + zlevel = pop!(attributes, :zlevel, 0) + # dest = pop!(attributes, :dest, nothing) plot = Plot{Makie.default_plot_func(F, args)}(args, attributes) - isnothing(source) || (plot.kw[:source] = source) - isnothing(dest) || (plot.kw[:dest] = dest) + plot.kw[:source] = source + plot.kw[:zlevel] = zlevel + # isnothing(dest) || (plot.kw[:dest] = dest) Makie.plot!(ax, plot) return plot end diff --git a/test/globeaxis.jl b/test/globeaxis.jl new file mode 100644 index 00000000..33fe035f --- /dev/null +++ b/test/globeaxis.jl @@ -0,0 +1,2 @@ +fig = Figure() +ga = GlobeAxis(fig[1,1]) diff --git a/test/globetransforms.jl b/test/globetransforms.jl new file mode 100644 index 00000000..e69de29b From a7aca84f565d23069ba26a414a5ee4c55f7f35f3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 20 Nov 2024 16:33:27 -0500 Subject: [PATCH 05/12] fix GeodesyGlobeTransform constructors --- src/sphere/globetransform.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/sphere/globetransform.jl b/src/sphere/globetransform.jl index 5882a4a7..9d51ae6f 100644 --- a/src/sphere/globetransform.jl +++ b/src/sphere/globetransform.jl @@ -57,16 +57,18 @@ struct GeodesyGlobeTransform{CoordType, T} <: GlobeTransform end -function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.LLA}) - return GeodesyGlobeTransform{src}(Geodesy.ECEFfromLLA(dest), zlevel) +function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.LLA}, zlevel) + transf = Geodesy.ECEFfromLLA(dest) + return GeodesyGlobeTransform{src, typeof(transf)}(transf, zlevel) end -function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.UTM}) - return GeodesyGlobeTransform{src}(Geodesy.ECEFfromUTM(dest), zlevel) +function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.UTM}, zlevel) + transf = Geodesy.ECEFfromUTM(dest) + return GeodesyGlobeTransform{src, typeof(transf)}(transf, zlevel) end -function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.UTMZ}) +function GeodesyGlobeTransform(dest::Geodesy.Ellipsoid, src::Type{Geodesy.UTMZ}, zlevel) return GeodesyGlobeTransform{src}(Geodesy.ECEFfromUTMZ(dest), zlevel) end From c4f61d3ca016d6d42796e13b00f6a6476b32622d Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 20 Nov 2024 18:31:41 -0500 Subject: [PATCH 06/12] apply 3d triangulation to GlobeTransform too --- src/triangulation3d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/triangulation3d.jl b/src/triangulation3d.jl index a6f914cb..d34fd7ed 100644 --- a/src/triangulation3d.jl +++ b/src/triangulation3d.jl @@ -26,7 +26,7 @@ end # We don't want to override the general case poly_convert for all polygons, because that's piracy, # but we _can_ override it for the specific case of a 3D polygon that is being transformed by a function # that is a subtype of Union{<: Proj.Transformation, <: GeoMakie.Geodesy.ECEFfromLLA} -function Makie.poly_convert(polygon::GeometryBasics.Polygon, transform_func::Union{<: Proj.Transformation, <: GeoMakie.Geodesy.ECEFfromLLA}) +function Makie.poly_convert(polygon::GeometryBasics.Polygon, transform_func::Union{<: Proj.Transformation, <: GeoMakie.Geodesy.ECEFfromLLA, <: GlobeTransform}) outer = GeometryBasics.metafree(GeometryBasics.coordinates(polygon.exterior)) PT = Makie.float_type(outer) From c9e5c3f2dfa76253beb3375aaca591ea77af4021 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 27 Nov 2024 14:11:36 -0500 Subject: [PATCH 07/12] add that to GlobeAxis too --- src/sphere/globeaxis.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/sphere/globeaxis.jl b/src/sphere/globeaxis.jl index 139b88a9..3d765af1 100644 --- a/src/sphere/globeaxis.jl +++ b/src/sphere/globeaxis.jl @@ -332,6 +332,7 @@ tightlimits!(::GlobeAxis) = nothing # TODO implement!? By getting the bbox of t # This is where we override the stuff to make it our stuff. function Makie.plot!(axis::GlobeAxis, plot::Makie.AbstractPlot) + # deal with setting the transform_func correctly source = pop!(plot.kw, :source, axis.source) zlevel = pop!(plot.kw, :zlevel, 0) # @show plot.kw @@ -339,6 +340,8 @@ function Makie.plot!(axis::GlobeAxis, plot::Makie.AbstractPlot) trans = Makie.Transformation(transformfunc; get(plot.kw, :transformation, Attributes())...) plot.kw[:transformation] = trans + reset_limits = to_value(pop!(plot.kw, :reset_limits, true)) + Makie.plot!(axis.scene, plot) # reassign popped observables back to attributes, so that they can be called by the user... @@ -348,7 +351,7 @@ function Makie.plot!(axis::GlobeAxis, plot::Makie.AbstractPlot) # adjust the limit margins in those cases automatically. # However, for spheres, we want to keep user zoom level if possible. # Makie.needs_tight_limits(plot) && Makie.tightlimits!(axis) - if Makie.is_open_or_any_parent(axis.scene) + if Makie.is_open_or_any_parent(axis.scene) && reset_limits Makie.reset_limits!(axis) end return plot From 8b6c6a648a36e9cfe328766d608f7db89915ce8d Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 23 Dec 2024 09:59:06 +0530 Subject: [PATCH 08/12] support GB 0.5 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c1f27545..a0f8f1d7 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,7 @@ GeoInterface = "0.5, 1.0" GeoInterfaceMakie = "0.1.6" GeoJSON = "0.6, 0.7, 0.8" Geodesy = "1.1.0" -GeometryBasics = "0.4.11" +GeometryBasics = "0.4.11, 0.5" GeometryOps = "0.1.6" ImageIO = "0.6" LinearAlgebra = "1" From 6eacfac636c927aeeeacaaa611ffdd12da87dc86 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 23 Dec 2024 09:59:18 +0530 Subject: [PATCH 09/12] pass all attributes through to mesh using shared_attributes --- src/mesh_image.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/mesh_image.jl b/src/mesh_image.jl index d581bd2e..496ad99a 100644 --- a/src/mesh_image.jl +++ b/src/mesh_image.jl @@ -114,23 +114,31 @@ function Makie.plot!(plot::MeshImage) # Let's remedy that now! final_mesh = lift(plot, points_observable, faces_observable, uv_observable; ignore_equal_values = true#=, priority = -100=#) do points, faces, uv return GeometryBasics.Mesh( - GeometryBasics.meta(points; uv=uv), # each point gets a UV, they're interpolated on faces - faces + points, + faces; + uv = uv, # each point gets a UV, they're interpolated on faces ) end # Finally, we plot the mesh. + shared_attrs = Makie.shared_attributes(plot, Makie.Mesh) + + haskey(shared_attrs, :color) && pop!(shared_attrs, :color) + haskey(shared_attrs, :shading) && pop!(shared_attrs, :shading) + haskey(shared_attrs, :transformation) && pop!(shared_attrs, :transformation) + haskey(shared_attrs, :uv_transform) && pop!(shared_attrs, :uv_transform) + mesh!( plot, final_mesh; color = plot.converted[3], # pass on the color directly - MakieCore.colormap_attributes(plot)..., # pass on all colormap attributes shading = NoShading, # transformation = Makie.Transformation( plot.transformation; # connect up the model matrix to the parent's model matrix transform_func = identity # do NOT connect the transform func, since we've already done that. identity provides a custom transform func, while `nothing` signals that you don't care. ), - uv_transform = plot.uv_transform + uv_transform = plot.uv_transform, + shared_attrs... ) # TODO: get a `:transformed` space out so we don't need this `transformation` hack end From 762bb62d41459c7918111279cadf0cdc436372c2 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 25 Dec 2024 12:45:44 +0530 Subject: [PATCH 10/12] Handle both versions of GB in meshimage recipe --- src/mesh_image.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/mesh_image.jl b/src/mesh_image.jl index 496ad99a..37c273e6 100644 --- a/src/mesh_image.jl +++ b/src/mesh_image.jl @@ -113,11 +113,18 @@ function Makie.plot!(plot::MeshImage) # You may have noticed that nowhere above did we actually create a mesh. # Let's remedy that now! final_mesh = lift(plot, points_observable, faces_observable, uv_observable; ignore_equal_values = true#=, priority = -100=#) do points, faces, uv - return GeometryBasics.Mesh( - points, - faces; - uv = uv, # each point gets a UV, they're interpolated on faces - ) + @static if hasproperty(GeometryBasics, :metafree) # pre v0.5 + return GeometryBasics.Mesh( + GeometryBasics.meta(points; uv=uv), # each point gets a UV, they're interpolated on faces + faces + ) + else # post v0.5 / Makie v0.22 + return GeometryBasics.Mesh( + points, + faces; + uv = uv, # each point gets a UV, they're interpolated on faces + ) + end end # Finally, we plot the mesh. From d31193eed0b7c89e7642501af0313b12111e34b4 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 26 Dec 2024 12:26:24 +0530 Subject: [PATCH 11/12] add Makie master versions --- .github/workflows/ci.yml | 7 ++++--- Project.toml | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index dc3f46b7..81b0350c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,9 +70,10 @@ jobs: Pkg.add([ PackageSpec(path = pwd()), PackageSpec(name = "DocumenterVitepress", rev = "master"), - # PackageSpec(name = "Makie", rev = "master"), - # PackageSpec(name = "MakieCore", rev = "master"), - # PackageSpec(name = "CairoMakie", rev = "master"), + PackageSpec(name = "Makie", rev = "master"), + PackageSpec(name = "MakieCore", rev = "master"), + PackageSpec(name = "CairoMakie", rev = "master"), + PackageSpec("GLMakie", rev = "master"), PackageSpec(name = "Rasters", rev = "main"), PackageSpec(name = "ArchGDAL", rev = "master"), PackageSpec(url = "https://github.com/asinghvi17/OhMyCards.jl", rev = "main"), diff --git a/Project.toml b/Project.toml index a0f8f1d7..ba3c7f00 100644 --- a/Project.toml +++ b/Project.toml @@ -36,7 +36,7 @@ GeometryBasics = "0.4.11, 0.5" GeometryOps = "0.1.6" ImageIO = "0.6" LinearAlgebra = "1" -Makie = "0.21.6" +Makie = "0.21.6, 0.22" NaturalEarth = "0.1" Proj = "1" Reexport = "1" From 817b2463fa2983e2e8546e996eccdb3bf2d87563 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 26 Dec 2024 12:26:47 +0530 Subject: [PATCH 12/12] fix typo --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 81b0350c..9d07ecb5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -73,7 +73,7 @@ jobs: PackageSpec(name = "Makie", rev = "master"), PackageSpec(name = "MakieCore", rev = "master"), PackageSpec(name = "CairoMakie", rev = "master"), - PackageSpec("GLMakie", rev = "master"), + PackageSpec(name = "GLMakie", rev = "master"), PackageSpec(name = "Rasters", rev = "main"), PackageSpec(name = "ArchGDAL", rev = "master"), PackageSpec(url = "https://github.com/asinghvi17/OhMyCards.jl", rev = "main"),