|
|
|
@ -1,5 +1,5 @@
|
|
|
|
#!/usr/bin/env julia
|
|
|
|
#!/usr/bin/env julia
|
|
|
|
#
|
|
|
|
|
|
|
|
using DataFrames
|
|
|
|
using DataFrames
|
|
|
|
using Dates
|
|
|
|
using Dates
|
|
|
|
using CSV
|
|
|
|
using CSV
|
|
|
|
@ -12,140 +12,157 @@ using Printf
|
|
|
|
using TimeZones
|
|
|
|
using TimeZones
|
|
|
|
|
|
|
|
|
|
|
|
function get_station_lat_lon(station_id)
|
|
|
|
function get_station_lat_lon(station_id)
|
|
|
|
stations = readxml("activestations.xml")
|
|
|
|
stations = readxml("activestations.xml")
|
|
|
|
station_node = findfirst("/stations/station[@id=\"$station_id\"]", stations)
|
|
|
|
station_node = findfirst("/stations/station[@id=\"$station_id\"]", stations)
|
|
|
|
return parse.(Float32, [station_node["lat"], station_node["lon"]])
|
|
|
|
return parse.(Float32, [station_node["lat"], station_node["lon"]])
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
tryparsem(T, str) = something(tryparse(T, str), missing)
|
|
|
|
tryparsem(T, str) = something(tryparse(T, str), missing)
|
|
|
|
|
|
|
|
|
|
|
|
function get_tile_xy_fraction(lat, lon, zoom)
|
|
|
|
function get_tile_xy_fraction(lat, lon, zoom)
|
|
|
|
n = 2^zoom
|
|
|
|
n = 2^zoom
|
|
|
|
lat_rad = lat * pi / 180
|
|
|
|
lat_rad = lat * pi / 180
|
|
|
|
xtile = (lon + 180) / 360 * n
|
|
|
|
xtile = (lon + 180) / 360 * n
|
|
|
|
ytile = (1 - log(tan(lat_rad) + 1 / cos(lat_rad)) / pi) / 2 * n
|
|
|
|
ytile = (1 - log(tan(lat_rad) + 1 / cos(lat_rad)) / pi) / 2 * n
|
|
|
|
return xtile, ytile
|
|
|
|
return xtile, ytile
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_tile_xy(lat, lon, zoom)
|
|
|
|
function get_tile_xy(lat, lon, zoom)
|
|
|
|
xf, yf = get_tile_xy_fraction(lat, lon, zoom)
|
|
|
|
xf, yf = get_tile_xy_fraction(lat, lon, zoom)
|
|
|
|
return (Int(floor(xf)), Int(floor(yf)))
|
|
|
|
return (Int(floor(xf)), Int(floor(yf)))
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_tile(x, y, zoom)
|
|
|
|
function get_tile(x, y, zoom)
|
|
|
|
url = "https://tile.openstreetmap.org/$zoom/$x/$y.png"
|
|
|
|
url = "https://tile.openstreetmap.org/$zoom/$x/$y.png"
|
|
|
|
tmp_path = download(url)
|
|
|
|
tmp_path = download(url)
|
|
|
|
tile = load(tmp_path)
|
|
|
|
tile = load(tmp_path)
|
|
|
|
return tile
|
|
|
|
return tile
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_tile_by_lat_lon(lat, lon, zoom)
|
|
|
|
function get_tile_by_lat_lon(lat, lon, zoom)
|
|
|
|
xtile, ytile = get_tile_xy(lat, lon, zoom)
|
|
|
|
xtile, ytile = get_tile_xy(lat, lon, zoom)
|
|
|
|
return get_tile(xtile, ytile, zoom)
|
|
|
|
return get_tile(xtile, ytile, zoom)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_meters_per_pixel(lat, zoom)
|
|
|
|
function get_meters_per_pixel(lat, zoom)
|
|
|
|
px_per_tile = 256
|
|
|
|
px_per_tile = 256
|
|
|
|
return 156543.03 * cos(lat) / (2^zoom)
|
|
|
|
return 156543.03 * cos(lat) / (2^zoom)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_map(lat, lon, zoom)
|
|
|
|
function get_map(lat, lon, zoom)
|
|
|
|
px_per_tile = 256
|
|
|
|
px_per_tile = 256
|
|
|
|
meters_per_pixel = get_meters_per_pixel(lat, zoom)
|
|
|
|
meters_per_pixel = get_meters_per_pixel(lat, zoom)
|
|
|
|
xf, yf = get_tile_xy_fraction(lat, lon, zoom)
|
|
|
|
xf, yf = get_tile_xy_fraction(lat, lon, zoom)
|
|
|
|
x0, y0 = (Int(floor(xf)), Int(floor(yf)))
|
|
|
|
x0, y0 = (Int(floor(xf)), Int(floor(yf)))
|
|
|
|
xpoint = 1 + px_per_tile + Int(floor((xf - x0) * px_per_tile))
|
|
|
|
xpoint = 1 + px_per_tile + Int(floor((xf - x0) * px_per_tile))
|
|
|
|
ypoint = 1 + px_per_tile + Int(floor((yf - y0) * px_per_tile))
|
|
|
|
ypoint = 1 + px_per_tile + Int(floor((yf - y0) * px_per_tile))
|
|
|
|
map = Array{RGB{N0f8}, 2}(undef, (px_per_tile * 3, px_per_tile * 3))
|
|
|
|
map = Array{RGB{N0f8}, 2}(undef, (px_per_tile * 3, px_per_tile * 3))
|
|
|
|
for i in -1:1
|
|
|
|
for i in -1:1
|
|
|
|
x = x0 + i
|
|
|
|
x = x0 + i
|
|
|
|
ioffset = 1 + (i + 1) * px_per_tile
|
|
|
|
ioffset = 1 + (i + 1) * px_per_tile
|
|
|
|
for j in -1:1
|
|
|
|
for j in -1:1
|
|
|
|
y = y0 + j
|
|
|
|
y = y0 + j
|
|
|
|
joffset = 1 + (j + 1) * px_per_tile
|
|
|
|
joffset = 1 + (j + 1) * px_per_tile
|
|
|
|
img = get_tile(x, y, zoom)
|
|
|
|
img = get_tile(x, y, zoom)
|
|
|
|
map[joffset:joffset+px_per_tile-1, ioffset:ioffset+px_per_tile-1] = img
|
|
|
|
map[joffset:joffset+px_per_tile-1, ioffset:ioffset+px_per_tile-1] = img
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)),
|
|
|
|
draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)),
|
|
|
|
RGB{N0f8}(0.5, 0.0, 0.5))
|
|
|
|
RGB{N0f8}(0.5, 0.0, 0.5))
|
|
|
|
return map, xpoint, ypoint
|
|
|
|
return map, xpoint, ypoint
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_realtime_data_file(station)
|
|
|
|
function get_realtime_data_file(station)
|
|
|
|
station = uppercase(station)
|
|
|
|
station = uppercase(station)
|
|
|
|
url = "http://www.ndbc.noaa.gov/data/realtime2/$station.txt"
|
|
|
|
url = "http://www.ndbc.noaa.gov/data/realtime2/$station.txt"
|
|
|
|
return download(url)
|
|
|
|
return download(url)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
|
|
function get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
|
|
station_fname = get_realtime_data_file(station_name)
|
|
|
|
station_fname = get_realtime_data_file(station_name)
|
|
|
|
lines = open(station_fname) do file
|
|
|
|
lines = open(station_fname) do file
|
|
|
|
readlines(file)
|
|
|
|
readlines(file)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
dt_array = (y, mon, d, hr, min) -> (astimezone.(ZonedDateTime.(
|
|
|
|
dt_array = (y, mon, d, hr, min) -> (astimezone.(ZonedDateTime.(
|
|
|
|
DateTime.(y, mon, d, hr, min), tz"UTC"), timezone))
|
|
|
|
DateTime.(y, mon, d, hr, min), tz"UTC"), timezone))
|
|
|
|
deg_c_to_f = passmissing(x -> x * 9.0 / 5.0 + 32.0)
|
|
|
|
deg_c_to_f = passmissing(x -> x * 9.0 / 5.0 + 32.0)
|
|
|
|
speed_ms_to_mph = passmissing(x -> 0.44704 * x)
|
|
|
|
speed_ms_to_mph = passmissing(x -> 0.44704 * x)
|
|
|
|
lines[1] = lstrip(lines[1], '#')
|
|
|
|
lines[1] = lstrip(lines[1], '#')
|
|
|
|
data = join([lines[1:1]; lines[3:end]], "\n")
|
|
|
|
data = join([lines[1:1]; lines[3:end]], "\n")
|
|
|
|
df = CSV.read(IOBuffer(data), DataFrame; header=1, delim=" ",
|
|
|
|
df = CSV.read(IOBuffer(data), DataFrame; header=1, delim=" ",
|
|
|
|
ignorerepeated=true)
|
|
|
|
ignorerepeated=true)
|
|
|
|
select!(df,
|
|
|
|
select!(df,
|
|
|
|
[:YY, :MM, :DD, :hh, :mm] => dt_array => :TIME,
|
|
|
|
[:YY, :MM, :DD, :hh, :mm] => dt_array => :TIME,
|
|
|
|
:WDIR => x -> tryparsem.(Int, x),
|
|
|
|
:WDIR => x -> tryparsem.(Int, x),
|
|
|
|
:WSPD => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
|
|
|
|
:WSPD => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
|
|
|
|
:GST => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
|
|
|
|
:GST => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
|
|
|
|
:PRES => x -> tryparsem.(Float32, x),
|
|
|
|
:PRES => x -> tryparsem.(Float32, x),
|
|
|
|
:ATMP => x -> deg_c_to_f.(tryparsem.(Float32, x)),
|
|
|
|
:ATMP => x -> deg_c_to_f.(tryparsem.(Float32, x)),
|
|
|
|
:WTMP => x -> deg_c_to_f.(tryparsem.(Float32, x));
|
|
|
|
:WTMP => x -> deg_c_to_f.(tryparsem.(Float32, x));
|
|
|
|
renamecols=false)
|
|
|
|
renamecols=false)
|
|
|
|
return df
|
|
|
|
return df
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function plot_station_data(station_name, station_map, station_x, station_y,
|
|
|
|
function plot_station_data(station_name, station_map, station_x, station_y,
|
|
|
|
df, i)
|
|
|
|
df, i)
|
|
|
|
fig = Figure(resolution=size(station_map))
|
|
|
|
fig = Figure(resolution=size(station_map))
|
|
|
|
ax = CairoMakie.Axis(fig[1, 1])
|
|
|
|
ax = CairoMakie.Axis(fig[1, 1])
|
|
|
|
hidespines!(ax)
|
|
|
|
hidespines!(ax)
|
|
|
|
hidedecorations!(ax)
|
|
|
|
hidedecorations!(ax)
|
|
|
|
plot_station_data_to_axis(ax, station_name, station_map,
|
|
|
|
plot_station_data_to_axis(ax, station_name, station_map,
|
|
|
|
station_x, station_y, df, i)
|
|
|
|
station_x, station_y, df, i)
|
|
|
|
return fig
|
|
|
|
return fig
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function plot_station_data_to_axis(ax, station_name, station_map, station_x, station_y,
|
|
|
|
function plot_station_data_to_axis(ax, station_name, station_map, station_x, station_y,
|
|
|
|
df, i)
|
|
|
|
df, i)
|
|
|
|
obs_time = df.TIME[i]
|
|
|
|
obs_time = df.TIME[i]
|
|
|
|
obs_time_s = Dates.format(obs_time, "e, u d H:MM")
|
|
|
|
obs_time_s = Dates.format(obs_time, "e, u d H:MM")
|
|
|
|
image!(ax, rotr90(station_map))
|
|
|
|
image!(ax, rotr90(station_map))
|
|
|
|
ax.title = "$station_name $obs_time_s"
|
|
|
|
ax.title = "$station_name $obs_time_s"
|
|
|
|
wind_rad = deg2rad(90 - df.WDIR[1])
|
|
|
|
wind_rad = deg2rad(90 - df.WDIR[1])
|
|
|
|
wind_mph = df.WSPD[i]
|
|
|
|
wind_mph = df.WSPD[i]
|
|
|
|
txt = @sprintf "Air: %.1f F, Water: %.1f F, Wind: %.1f mph" df.ATMP[i] df.WTMP[i] wind_mph
|
|
|
|
if ismissing(wind_mph) || ismissing(wind_rad)
|
|
|
|
wind_vec = 50 * Vec2f(cos(wind_rad), sin(wind_rad)) * wind_mph / 10
|
|
|
|
wind_mph = 0.0
|
|
|
|
tooltip!(ax, station_x, station_y, txt; offset=10)
|
|
|
|
wind_rad = 0
|
|
|
|
arrows!(ax, [Point2f(station_x, station_y)], [wind_vec], fxaa=true,
|
|
|
|
end
|
|
|
|
linewidth=2, align=:center, arrowhead='⟰')
|
|
|
|
atmp = df.ATMP[i]
|
|
|
|
return ax
|
|
|
|
wtmp = df.WTMP[i]
|
|
|
|
|
|
|
|
if ismissing(atmp)
|
|
|
|
|
|
|
|
atmp = 0.0
|
|
|
|
|
|
|
|
end
|
|
|
|
|
|
|
|
if ismissing(wtmp)
|
|
|
|
|
|
|
|
wtmp = 0.0
|
|
|
|
|
|
|
|
end
|
|
|
|
|
|
|
|
txt = @sprintf "Air: %.1f F, Water: %.1f F, Wind: %.1f mph" atmp wtmp wind_mph
|
|
|
|
|
|
|
|
wind_vec = 50 * Vec2f(cos(wind_rad), sin(wind_rad)) * wind_mph / 10
|
|
|
|
|
|
|
|
tooltip!(ax, station_x, station_y, txt; offset=10)
|
|
|
|
|
|
|
|
arrows!(ax, [Point2f(station_x, station_y)], [wind_vec], fxaa=true,
|
|
|
|
|
|
|
|
linewidth=2, align=:center, arrowhead='⟰')
|
|
|
|
|
|
|
|
return ax
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
function save_animation(station_name, out_path, days, real_hour_per_animation_second=1)
|
|
|
|
function save_animation(station_name, out_path, hours, real_hour_per_animation_second)
|
|
|
|
lat, lon = get_station_lat_lon(station_name)
|
|
|
|
lat, lon = get_station_lat_lon(station_name)
|
|
|
|
station_map, station_x, station_y = get_map(lat, lon, 14)
|
|
|
|
station_map, station_x, station_y = get_map(lat, lon, 14)
|
|
|
|
df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
|
|
df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
|
|
diff_hours = (df.TIME[0] - df.TIME[1]).value / 1000 / 3600
|
|
|
|
diff_hours = (df.TIME[1] - df.TIME[2]).value / 1000 / 3600
|
|
|
|
fps = real_hour_per_animation_second / diff_hours
|
|
|
|
stride = Int(ceil(real_hour_per_animation_second / diff_hours))
|
|
|
|
nframes = fps * days * 24 * 3600
|
|
|
|
println("diff hours ", diff_hours)
|
|
|
|
|
|
|
|
println("data size ", size(df, 1))
|
|
|
|
|
|
|
|
println("stride ", stride)
|
|
|
|
|
|
|
|
end_index = Int(min(stride * hours + 1, size(df, 1)))
|
|
|
|
|
|
|
|
data_indexes = range(end_index, 1; step=-stride)
|
|
|
|
|
|
|
|
println("data indexes ", data_indexes)
|
|
|
|
|
|
|
|
|
|
|
|
fig = Figure(resolution=size(station_map))
|
|
|
|
fig = Figure(resolution=size(station_map))
|
|
|
|
ax = Axis(fig[1, 1])
|
|
|
|
ax = CairoMakie.Axis(fig[1, 1])
|
|
|
|
hidespines!(ax)
|
|
|
|
hidespines!(ax)
|
|
|
|
hidedecorations!(ax)
|
|
|
|
hidedecorations!(ax)
|
|
|
|
|
|
|
|
|
|
|
|
record(fig, out_path, range(1, nframes); framerate=fps) do i
|
|
|
|
record(fig, out_path, data_indexes; framerate=1) do i
|
|
|
|
plot_station_data_to_axis(ax, station_name, station_map,
|
|
|
|
plot_station_data_to_axis(ax, station_name, station_map,
|
|
|
|
station_x, station_y, df, nframes - i + 1)
|
|
|
|
station_x, station_y, df, end_index - i + 1)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
@ -157,16 +174,20 @@ function main()
|
|
|
|
station_name = ARGS[1]
|
|
|
|
station_name = ARGS[1]
|
|
|
|
|
|
|
|
|
|
|
|
out_fname = "$station_name.png"
|
|
|
|
out_fname = "$station_name.png"
|
|
|
|
if (length(ARGS) > 1)
|
|
|
|
if length(ARGS) > 1
|
|
|
|
out_fname = ARGS[2]
|
|
|
|
out_fname = ARGS[2]
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
lat, lon = get_station_lat_lon(station_name)
|
|
|
|
if any(endswith.(out_fname, [".gif", ".webm", ".mp4", ".mkv"]))
|
|
|
|
station_map, station_x, station_y = get_map(lat, lon, 14)
|
|
|
|
save_animation(station_name, out_fname, 24, 1)
|
|
|
|
df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
|
|
else
|
|
|
|
fig = plot_station_data(station_name, station_map, station_x, station_y,
|
|
|
|
lat, lon = get_station_lat_lon(station_name)
|
|
|
|
df, 1)
|
|
|
|
station_map, station_x, station_y = get_map(lat, lon, 14)
|
|
|
|
save("$station_name.png", fig)
|
|
|
|
df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
|
|
|
|
|
|
fig = plot_station_data(station_name, station_map, station_x, station_y,
|
|
|
|
|
|
|
|
df, 1)
|
|
|
|
|
|
|
|
save(out_fname, fig)
|
|
|
|
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
main()
|
|
|
|
main()
|
|
|
|
|