diff --git a/docs/Project.toml b/docs/Project.toml index 9058553e..060aa6fd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +Animations = "27a7e980-b3e6-11e9-2bcd-0b925532e340" ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" diff --git a/docs/make.jl b/docs/make.jl index 940c20f0..fe712dcf 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -37,6 +37,7 @@ makedocs(; "Basic examples" => "examples/basic.md", "New API" => "examples/new.md", "Orthographic projection" => "examples/orthographic.md", + "German Lakes" => "examples/german_lakes.md", "Geostationary satellite image" => "examples/geostationary_image.md", "Contourf" => "examples/contourf.md", # "Multiple CRS" => "examples/multiple_crs.md", diff --git a/examples/german_lakes.jl b/examples/german_lakes.jl new file mode 100644 index 00000000..6a4c336b --- /dev/null +++ b/examples/german_lakes.jl @@ -0,0 +1,14 @@ +# # German lakes +using GeoMakie, CairoMakie +using GeoJSON, Downloads +CairoMakie.activate!(px_per_unit = 4) # hide + +geoger = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/2_bundeslaender/4_niedrig.geo.json"), String)) +lakes = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_10m_lakes_europe.geojson"), String)) + +fig = Figure() +ga = GeoAxis(fig[1, 1]; dest = "+proj=merc", limits=((6, 15), (47, 55))) + +poly!(ga, geoger; strokewidth = 0.7, color=:gold, rasterize = 5) +poly!(ga, lakes; strokewidth = 0.7, color=:blue, rasterize = 5, xautolimits=false, yautolimits=false) +fig diff --git a/examples/is_it_a_plane.jl b/examples/is_it_a_plane.jl new file mode 100644 index 00000000..f665962c --- /dev/null +++ b/examples/is_it_a_plane.jl @@ -0,0 +1,59 @@ +# # Geodesic paths - animation + +# Let's take the great circle flight path from New York (JFK) +# to Singapore (SIN) airport. + +using GeoMakie, CairoMakie +CairoMakie.activate!(px_per_unit = 2) # hide +using Proj, Animations + +jfk = Point2f(-73.7789, 40.6397) +sin = Point2f(103.9894, 1.3592) + +# First, we define the globe, as the WGS84 ellipsoid: +geod = Proj.geod_geodesic(6378137, 1/298.257223563) +# Then, we can solve the inverse geodesic problem, which provides +# the shortest path between two points on our defined ellipsoid: +inv_line = Proj.geod_inverseline(geod, reverse(jfk)..., reverse(sin)...) +# Just for reference, this is the path: +f, a, p = lines(reverse(Proj.geod_path(geod, reverse(jfk)..., reverse(sin)...))...; linewidth = 3, axis = (; type = GeoAxis, dest = "+proj=natearth")); lines!(a, GeoMakie.coastlines(), color = (:black, 0.4)); f + +# We'll use a satellite view for this, and alter the projection as a way of controlling the animation. + +# First, we'll create 2 observables which control the position of the "camera": +# distance along path (from 0 to 1) and altitude (in meters)! + +# The projection will always be centered at wherever the plane is. + +# We first + +times = [0, 0.5, 17.5, 18] +distances = [0, 0.05, 0.95, 1] +altitudes = [357860, 35786000/2, 35786000/2, 357860] +distance_animation = Animation(times, distances, linear()) +altitude_animation = Animation(times, altitudes, sineio()) + + +# In order to investigate this kind of projection, you can create +# a GeoAxis with the projection you want, and then change the +# altitude to see how the zoom works in real time! +@time begin +fig = Figure() +sl = Slider(fig[2, 1], range = exp.(LinRange(log(357860), log(35786000), 30)), startvalue = 35786000) +satview_projection = lift(sl.value) do alt + "+proj=geos +h=$(round(Int, alt)) +lon_0=$(sin[1]) +lat_0=$(sin[2])" +end +ga = GeoAxis(fig[1, 1]; dest = satview_projection) +meshimage!(ga, -180..180, -90..90, GeoMakie.earth(), shading = NoShading) +fig +end + +record(fig, "plots/plane.mp4", LinRange(0, 1, 120)) do i + satview_projection[] = "+proj=geos +h=$(round(Int, at(altitude_animation, i*18))) +lon_0=$(Proj.geod_position_relative(inv_line, i)[2]) +lat_0=$(Proj.geod_position_relative(inv_line, i)[1])" + yield() +end + +fig = Figure() +ga = GeoAxis(fig[1, 1]; dest = "+proj=nsper +h=3000000 +lat_0=-20 +lon_0=145") +meshimage!(ga, -180..180, -90..90, GeoMakie.earth(), shading = NoShading) +fig \ No newline at end of file