-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtec_map_2D.jl
139 lines (122 loc) · 5.26 KB
/
tec_map_2D.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
using TEC_viz
using Dates
using GeoMakie
using GLMakie
using HDF5
## Extract the data
filename = joinpath("data", "gps150317g.004.hdf5")
# filename = joinpath(@__DIR__, "data", "gps131107g.004.hdf5")
fid = h5open(filename, "r")
data = read(fid)
close(fid)
lat = data["Data"]["Array Layout"]["gdlat"]
lon = data["Data"]["Array Layout"]["glon"]
timestamps = data["Data"]["Array Layout"]["timestamps"]
tec = data["Data"]["Array Layout"]["2D Parameters"]["tec"]
dtec = data["Data"]["Array Layout"]["2D Parameters"]["dtec"]
## Plot
# Prepare the figure and axis
fig = Figure(; size = (1500, 900))
ax1 = GeoAxis(fig[1, 1]; dest = "+proj=ortho +lat_0=90", xticks = 0:60:360,
yticks = 0:15:90,
xticklabelsvisible = false, yticklabelsvisible = false)
ax2 = GeoAxis(fig[1, 3]; dest = "+proj=ortho +lat_0=-90", xticks = 0:60:360,
yticks = 0:-15:-90,
xticklabelsvisible = false, yticklabelsvisible = false)
# Plot the coastlines and terrain
l1 = lines!(ax1, GeoMakie.coastlines(); color = :black)
l2 = lines!(ax2, GeoMakie.coastlines(); color = :black)
translate!(l1, 0, 0, 1) # put the coastlines over everything
translate!(l2, 0, 0, 1) # put the coastlines over everything
earth1 = surface!(ax1, -180 .. 180, 0 .. 90, zeros(180, 360);
color = rotr90(GeoMakie.earth()[1:180, :]))
earth2 = surface!(ax2, -180 .. 180, -90 .. 0, zeros(180, 360);
color = rotr90(GeoMakie.earth()[181:end, :]))
translate!(earth1, 0, 0, -1) # put the terrain under everything
translate!(earth2, 0, 0, -1) # put the terrain under everything
# Initialize Observables
tec_points = tec[1, :, :]
good_idx = findall(!isnan, tec_points)
good_lon = Observable([lon[i[1]] for i in good_idx])
good_lat = Observable([lat[i[2]] for i in good_idx])
good_tec = Observable(tec_points[good_idx])
# Plot the tec points
sc1 = scatter!(ax1, good_lon, good_lat; color = good_tec, colormap = :inferno,
colorrange = (0, 30))
# colorrange = (-1, 10))
sc2 = scatter!(ax2, good_lon, good_lat; color = good_tec, colormap = :inferno,
colorrange = (0, 30))
# colorrange = (-1, 10))
# Add a colobar
Colorbar(fig[1, 2], sc1; label = "TEC units", tellheight = false,
height = @lift Fixed($(pixelarea(ax1.scene)).widths[2]))
Colorbar(fig[1, 4], sc2; label = "TEC units", tellheight = false,
height = @lift Fixed($(pixelarea(ax1.scene)).widths[2]))
# Add a shade over the night side of Earth
nightshade_lon = Observable(night_shade(timestamps[1])[2])
surface!(ax1, nightshade_lon, 0 .. 90, ones(30, 15); color = fill((:black, 0.3), 30, 15))
surface!(ax2, nightshade_lon, -90 .. 0, ones(30, 15); color = fill((:black, 0.3), 30, 15))
# Rotate Earth so that the zenith is at the "top"
midnight_lon = Observable(night_shade(timestamps[1])[1])
ax1.dest[] = "+proj=ortho +lat_0=90 +lon_0=$(midnight_lon[])"
ax2.dest[] = "+proj=ortho +lat_0=-90 +lon_0=$(180+midnight_lon[])"
# Add a title
ax1.title = string(unix2datetime(timestamps[1]))
ax1.titlesize = 20
# Make a function to update the data points
function update_plot!(i_t, ax1, ax2, midnight_lon, nightshade_lon, tec, good_tec,
timestamps)
# update the Sun pointing direction and the nightshade accordingly
midnight_lon[], nightshade_lon[] = night_shade(timestamps[i_t])
ax1.dest[] = "+proj=ortho +lat_0=90 +lon_0=$(midnight_lon[])"
ax2.dest[] = "+proj=ortho +lat_0=-90 +lon_0=$(180+midnight_lon[])"
# update the tec points
local tec_points = tec[i_t, :, :]
local good_idx = findall(!isnan, tec_points)
good_lon.val = [lon[i[1]] for i in good_idx]
good_lat.val = [lat[i[2]] for i in good_idx]
good_tec[] = tec_points[good_idx]
notify(good_lon)
notify(good_lat)
# update the title
ax1.title = string(unix2datetime(timestamps[i_t]))
return nothing
end
# ## Testing if it works
# display(fig)
# for i_t in 1:length(timestamps)
# update_plot!(i_t, ax1, ax2, midnight_lon, nightshade_lon, tec, good_tec, timestamps)
# sleep(0.05)
# end
# Add a button to launch the animation
fig[2, :] = buttongrid = GridLayout(; tellwidth = false)
run = buttongrid[1, 1] = Button(fig; label = "Play/Pause", tellwidth = false)
isrunning = Observable(false)
i_t = Observable(1)
on(run.clicks) do clicks
isrunning[] = !isrunning[]
end
on(run.clicks) do clicks
@async while isrunning[]
i_t[] < length(timestamps) ? i_t[] += 1 : i_t[] = 1 # take next time step and loop when t_max is reached
isopen(fig.scene) || break # ensures animation stops if the figure is closed
update_plot!(i_t[], ax1, ax2, midnight_lon, nightshade_lon, tec, good_tec,
timestamps)
sleep(0.05)
end
end
# Add a button to reset the animation
reset = buttongrid[1, 2] = Button(fig; label = "Reset", tellwidth = false)
on(reset.clicks) do clicks
i_t[] = 1
update_plot!(i_t[], ax1, ax2, midnight_lon, nightshade_lon, tec, good_tec, timestamps)
end
# Finally display the figure
display(fig)
# display(GLMakie.Screen(), fig)
## That's to animate and save
display(fig)
video_file = joinpath(@__DIR__, "animations", "tec_map_20131107.mp4")
record(fig, video_file, 1:length(timestamps); px_per_unit = 2, framerate = 15) do i_t
update_plot!(i_t, ax1, ax2, midnight_lon, nightshade_lon, tec, good_tec, timestamps)
end