Skip to content

Commit

Permalink
Add German lakes example (#200)
Browse files Browse the repository at this point in the history
  • Loading branch information
asinghvi17 authored Feb 6, 2024
1 parent 972f7ac commit 445970c
Show file tree
Hide file tree
Showing 4 changed files with 75 additions and 0 deletions.
1 change: 1 addition & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
14 changes: 14 additions & 0 deletions examples/german_lakes.jl
Original file line number Diff line number Diff line change
@@ -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
59 changes: 59 additions & 0 deletions examples/is_it_a_plane.jl
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 445970c

Please sign in to comment.