Empty/deformed Polygons with LScene #282
Comments
|
I'm not sure what to do here... Can you make a minimal example that shows the problem? |
|
Yes, trying to boil this down, the most minimal I came up with so far is this using GLMakie, GeoMakie
import GeoMakie: GeometryBasics
import GeoMakie.Makie.GeometryBasics: Polygon
import LinearAlgebra: cross, dot
# Anshul defined this
function GeometryBasics.earcut_triangulate(polygon::Vector{Vector{Point{3,Float64}}})
# Here, we are assuming that the polygon is actually planar,
# but the plane is in 3D, and not necessarily the XY plane.
# So, we can find the plane of best fit using the first three points of the polygon!
p1, p2, p3 = polygon[1][1], polygon[1][2], polygon[1][3]
v1 = p2 - p1
v2 = p3 - p1
normal = cross(v1, v2)
x = v1
y = cross(normal, x)
projected_polygon = map(ring -> map(p -> Point2{Float64}(dot(p, x), dot(p, y)), ring), polygon)
lengths = map(x -> UInt32(length(x)), projected_polygon)
len = UInt32(length(lengths))
array = ccall((:u32_triangulate_f64, GeometryBasics.libearcut), Tuple{Ptr{GeometryBasics.GLTriangleFace},Cint},
(Ptr{Ptr{Float64}}, Ptr{UInt32}, UInt32), projected_polygon, lengths, len)
return unsafe_wrap(Vector{GeometryBasics.GLTriangleFace}, array[1], array[2]; own=true)
end
# coordinates of the vertices
londs = [45.0 135.0 225.0 315.0 0.0 90.0 90.0 180.0 180.0 270.0 270.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 45.0 135.0 225.0 315.0;
90.0 180.0 270.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 90.0 180.0 270.0 0.0;
45.0 135.0 225.0 315.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 90.0 90.0 180.0 180.0 270.0 270.0 0.0 45.0 135.0 225.0 315.0;
0.0 90.0 180.0 270.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 90.0 180.0 270.0;
45.0 135.0 225.0 315.0 0.0 90.0 90.0 180.0 180.0 270.0 270.0 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0 0.0 45.0 135.0 225.0 315.0]
latds = [90.0 90.0 90.0 90.0 59.4444 59.4444 59.4444 59.4444 59.4444 59.4444 59.4444 59.4444 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 -19.8757 -19.8757 -19.8757 -19.8757;
59.4444 59.4444 59.4444 59.4444 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -59.4444 -59.4444 -59.4444 -59.4444;
19.8757 19.8757 19.8757 19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -59.4444 -59.4444 -59.4444 -59.4444 -59.4444 -59.4444 -59.4444 -59.4444 -90.0 -90.0 -90.0 -90.0;
59.4444 59.4444 59.4444 59.4444 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -19.8757 -59.4444 -59.4444 -59.4444 -59.4444;
90.0 90.0 90.0 90.0 59.4444 59.4444 59.4444 59.4444 59.4444 59.4444 59.4444 59.4444 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 19.8757 -19.8757 -19.8757 -19.8757 -19.8757]
# create the polygons from the faces
faces = [Point2(λ, φ) for (λ, φ) in zip(londs, latds)]
polygons = [Polygon(faces[:, ij]) for ij in axes(faces, 2)]
transf = GeoMakie.Geodesy.ECEFfromLLA(GeoMakie.Geodesy.WGS84())
fig = Figure(size=(800, 800));
ax = LScene(fig[1,1], show_axis=false);
p = poly!(ax, polygons, color=rand(length(polygons)))
p.transformation.transform_func[] = transf
c = lines!(GeoMakie.coastlines(50); color=:black)
c.transformation.transform_func[] = transf
cc = cameracontrols(ax.scene)
cc.settings.mouse_translationspeed[] = 0.0
cc.settings.zoom_shift_lookat[] = false
Makie.update_cam!(ax.scene, cc)
fig
The bright yellow polygon in the south atlantic is not fully drawn for some reason I don't understand |
|
Even more minimal example, if you take polygon 12 and 20 of those 24 defined here they should have vertices located at julia> faces[:, 12]
5-element Vector{Point{2, Float64}}:
[0.0, 59.4444]
[0.0, 19.8757]
[315.0, -19.8757]
[315.0, 19.8757] # this one is missing
[0.0, 59.4444]
julia> faces[:, 20]
5-element Vector{Point{2, Float64}}:
[0.0, 19.8757] # this one
[0.0, -19.8757]
[0.0, -59.4444]
[315.0, -19.8757]
[0.0, 19.8757] # and this one are missing
|
|
Looks to me like the |
|
Yes, here you can see how the yellow diamond at the horizon in Russia is correctly split into a north and a south triangle but around the equator it seems to try to split into an east and west triangle which seems to "bury" half the polygon underneath the Earth's surface or similar?
|
|
|
I'm unintentionally creating some empty/deformed polygons with LScene that occur only around the equator for polygons that touch/cross the equator. The clear cut of the southern hemisphere is the equator but the empty polygons allow one to look through the sphere such that the same issue on the other side of the globe becomes visible
this happens in the
mk/globebranch in this pull request SpeedyWeather/SpeedyWeather.jl#600.the vertices themselves are fine
The text was updated successfully, but these errors were encountered: